#rm(list = ls(all = TRUE))
#install.packages("INLA", repos="http://www.math.ntnu.no/inla/R/testing")
#install.packages("bigmemory")
#install.packages("snow")
#install.packages("Rmpi")
#install.packages("ade4")
#install.packages("sp")
#install.packages("BAS")
#install.packages("hash")
#install.packages("stringi")
library(hash)
library(sp)
library(INLA)
library(parallel)
library(bigmemory)
library(snow)
library(MASS)
library(ade4)
#library(copula)
#library(compiler)
library(BAS)
library(stringi)
require(stats)
#compile INLA
estimate.bas.glm <- function(formula, data, family, prior, logn)
{
#only poisson and binomial families are currently adopted
X <- model.matrix(object = formula,data = data)
out <- bayesglm.fit(x = X, y = data[,1], family=family,coefprior=prior)
# use dic and aic as bic and aic correspondinly
return(list(mlik = out$logmarglik,waic = -(out$deviance + 2*out$rank) , dic = -(out$deviance + logn*out$rank)))
}
parallelize<-function(X,FUN)
{
max.cpu <- length(X)
cl <-makeCluster(max.cpu ,type = paral.type,outfile = "")#outfile = ""
clusterEvalQ(cl = cl,expr = c(library(INLA),library(bigmemory)))
if(exists("statistics"))
{
clusterExport(cl=cl, "statistics")
clusterEvalQ(cl,{statistics <- attach.resource(statistics);1})
}
res.par <- parLapply(cl = cl, X, FUN)
stopCluster(cl)
return(res.par)
}
estimate.bas.lm <- function(formula, data, prior, n, g = 0)
{
out <- lm(formula = formula,data = data)
# 1 for aic, 2 bic prior, else g.prior
p <- out$rank
if(prior == 1)
{
ss<-sum(out$residuals^2)
logmarglik <- -0.5*(log(ss)+2*p)
}
else if(prior ==2)
{
ss<-sum(out$residuals^2)
logmarglik <- -0.5*(log(ss)+log(n)*p)
}
else
{
Rsquare <- summary(out)$r.squared
#logmarglik = .5*(log(1.0 + g) * (n - p -1) - log(1.0 + g * (1.0 - Rsquare)) * (n - 1))*(p!=1)
logmarglik = .5*(log(1.0 + g) * (n - p) - log(1.0 + g * (1.0 - Rsquare)) * (n - 1))*(p!=1)
}
# use dic and aic as bic and aic correspondinly
return(list(mlik = logmarglik,waic = AIC(out) , dic = BIC(out)))
}
# covariates are thus
# 1) A constant +1
# 2) Base structures (CHG, CGH, CHH) +2
# 3) Gene class (5 classes) +4
# 4) Position within a gene (4 positions including promoters, after coding parts, intrones, exones) +3
# 5) Gene expression (continious) +1
# 6) A CpG island factor +1
# 7) A tranposone factor +1
# 8) combinations (if reasonable)
#closeAllConnections()
EMJMCMC2016 <- setRefClass(Class = "EMJMCMC2016",
fields = list(estimator.args = "list",
max.cpu = "integer",
objective = "integer",
p.prior = "numeric",
min.N = "integer",
estimator = "function",
parallelize = "function",
parallelize.global = "function",
parallelize.hyper = "function",
min.N.randomize = "integer",
max.N.randomize = "integer",
type.randomize = "integer",
thin_rate = "integer",
aa = "numeric",
cc = "numeric",
printable.opt = "logical",
fobserved = "vector",
switch.type = "integer",
n.size ="integer",
LocImprove = "array",
max.N = "integer",
recalc.margin = "numeric",
max.N.glob = "integer",
min.N.glob = "integer",
max.cpu.glob = "integer",
max.cpu.hyper = "integer",
switch.type.glob = "integer",
isobsbinary = "array",
fparam = "vector",
p.add = "array",
latent.formula = "character",
Nvars = "integer",
seed = "integer",
M.nd = "integer",
locstop.nd = "logical",
M.mcmc = "integer",
SA.param = "list",
p.allow.tree = "numeric",
p.nor = "numeric",
p.and = "numeric",
p.allow.replace = "numeric",
max.tree.size = "integer",
double.hashing = "logical",
hash.length = "integer",
update.marg.mc = "logical",
Nvars.max = "integer",
Nvars.init = "integer",
allow_offsprings = "integer",
mutation_rate = "integer",
g.results = "big.matrix"
),
methods = list(
#class constructor
initialize = function(estimator.function = inla, estimator.args.list = list(family = "gaussian",data = data.example, control.fixed=list(prec=list(default= 0.00001),prec.intercept = 0.00001,mean=list(default= 0),mean.intercept = 0),
control.family = list(hyper = list(prec = list(prior = "loggamma",param = c(0.00001,0.00001),initial = 0))),
control.compute = list(dic = TRUE, waic = TRUE, mlik = TRUE)), search.args.list = NULL,latent.formula = "")
{
estimator <<- estimator.function
estimator.args <<- estimator.args.list
latent.formula <<- latent.formula
g.results <<- big.matrix(nrow = 4,ncol = 2)
g.results[1,1]<- -Inf
g.results[1,2]<- 1
g.results[2,1]<- Inf
g.results[2,2]<- 1
g.results[3,1]<- Inf
g.results[3,2]<- 1
g.results[4,1]<- 0
g.results[4,2]<- 0
if(is.null(search.args.list))
{
max.cpu <<- as.integer(Nvars*0.05 + 1)
objective <<- as.integer(1)
if(Sys.info()['sysname']=="Windows")
{
parallelize <<- lapply
parallelize.global <<- lapply
parallelize.hyper <<- lapply
}
else
{
parallelize <<- mclapply
parallelize.global <<- mclapply
parallelize.hyper <<- mclapply
}
Nvars <<- as.integer(length(fparam.example))
min.N <<- as.integer(Nvars/6)
min.N.glob <<- as.integer(Nvars/3)
max.N.glob <<- as.integer(Nvars/2)
max.N <<- as.integer(Nvars/5)
switch.type.glob <<- as.integer(2)
min.N.randomize <<- as.integer(1)
max.N.randomize <<- as.integer(1)
type.randomize <<- as.integer(3)
max.cpu.glob <<- as.integer(Nvars*0.05 + 1)
max.cpu.hyper <<- as.integer(2)
printable.opt <<- FALSE
thin_rate<<- as.integer(-1)
p.allow.tree <<- 0.6
p.allow.replace <<- 0.3
max.tree.size<<- as.integer(15)
Nvars.max <<- as.integer(Nvars+10)
Nvars.init <<- as.integer(Nvars)
allow_offsprings <<- as.integer(0)
mutation_rate <<- as.integer(100)
locstop.nd <<-FALSE
double.hashing <<- (Nvars > 20)
hash.length <<- as.integer(25)
aa <<- 0.9
cc <<- 0.0
M.nd <<- as.integer(Nvars)
M.mcmc <<- as.integer(5)
SA.param <<- list(t.min = 0.0001,t.init = 10, dt = 3, M = as.integer(Nvars/5+1))
fobserved <<- fobserved.example
switch.type <<- as.integer(2)
n.size <<- as.integer(10)
LocImprove <<- as.array(c(50,50,50,50,150))
isobsbinary <<- as.array(0:(length(fparam.example)-1))
fparam <<- fparam.example
p.add <<- array(data = 0.5,dim = Nvars)
p.nor <<- 0.3
p.and <<- 0.7
if(exists("statistics"))
{
recalc.margin <<- 2^Nvars
}else if(exists("statistics1"))
{
recalc.margin <<- 2^Nvars
}else
{
recalc.margin <<- 2^Nvars
}
seed <<- as.integer(runif(n = 1,min = 1,max = 10000))
p.prior <<- runif(n = Nvars, min = 0.5,max = 0.5)
}
else
{
max.cpu <<- as.integer(search.args.list$max.cpu)
objective <<- as.integer(search.args.list$objective)
parallelize <<- search.args.list$parallelize
parallelize.global <<- search.args.list$parallelize.global
parallelize.hyper <<- search.args.list$parallelize.hyper
p.prior <<- search.args.list$p.prior
min.N <<- as.integer(search.args.list$min.N)
printable.opt <<- search.args.list$printable.opt
min.N.glob <<- as.integer(search.args.list$min.N.glob)
max.N.glob <<- as.integer(search.args.list$max.N.glob)
switch.type.glob <<- as.integer(search.args.list$switch.type.glob)
min.N.randomize <<- as.integer(search.args.list$min.N.randomize)
max.N.randomize <<- as.integer(search.args.list$max.N.randomize)
type.randomize <<- as.integer(search.args.list$type.randomize)
max.cpu.glob <<- as.integer(search.args.list$max.cpu.glob)
locstop.nd <<-search.args.list$locstop.nd
max.cpu.hyper <<- as.integer(search.args.list$max.cpu.hyper)
aa <<- search.args.list$lambda.a
thin_rate <-search.args.list$thin_rate
cc <<- search.args.list$lambda.c
M.nd <<- as.integer(search.args.list$stepsGreedy)
M.mcmc <<- as.integer(search.args.list$stepsLocMCMC)
SA.param <<- search.args.list$SA.params
fobserved <<- search.args.list$fobserved
switch.type <<- as.integer(search.args.list$fswitch.type)
n.size <<- as.integer(search.args.list$n.size)
LocImprove <<-search.args.list$prior.optimizer.freq
max.N <<- as.integer(search.args.list$max.N)
fparam <<- search.args.list$fparam
isobsbinary <<- as.array(0:(length(fparam)-1))
p.add <<- as.array(search.args.list$p.add)
recalc.margin <<- search.args.list$recalc.margin
Nvars <<- as.integer(length(fparam))
seed <<- search.args.list$seed
max.tree.size<<- as.integer(search.args.list$max.tree.size)
Nvars.max <<- as.integer(search.args.list$Nvars.max)
Nvars.init <<- as.integer(search.args.list$Nvars)
allow_offsprings <<- as.integer(search.args.list$allow_offsprings)
mutation_rate <<- as.integer(search.args.list$mutation_rate)
p.allow.tree <<- search.args.list$p.allow.tree
p.allow.replace <<- search.args.list$p.allow.replace
p.nor <<- search.args.list$p.nor
p.and <<- search.args.list$p.and
double.hashing <<- search.args.list$double.hashing
hash.length <<- as.integer(search.args.list$hash.length)
}
},
#transform binary numbers to decimal
bittodec.alt = function(bit) #transform a binary vector into a natural number to correspond between vector of solutions and storage array
{
n<-length(bit)
dec <- 0
for(i in 1:n)
{
j<-n-i
dec <- dec + ((2)^j)*bit[i]
}
return(dec)
},
bittodec = function(bit) #transform a binary vector into a natural number to correspond between vector of solutions and storage array
{
if(!double.hashing){
n<-length(bit)
dec <- 0
for(i in 1:n)
{
j<-n-i
dec <- dec + ((2)^j)*bit[i]
}
return(dec)
}
else
{
if(exists("statistics1"))
{
hash.level<-0
dec<- hashing(bit)+1
#print(dec)
sum.one<- sum(bit)*which.max(bit) + sum(bit[Nvars -7:Nvars])*Nvars
jjj<-1
while(!add.key(dec,bit,sum.one,T))
{
hash.level<-hash.level+1
bit1<- dectobit.alt(2654435761*(dec+97*sum.one+hash.level*36599)+hash.level*59+hash.level)
dec<- hashing(bit1)+1
#jjj<-jjj+1
#print(dec)
}
#print(jjj)
dec <- dec - 1
}
else if(exists("statistics"))
{
hash.level<-0
dec<- hashing(bit)+1
sum.one<- sum(bit)*which.max(bit) + sum(bit[Nvars -4:Nvars])*Nvars
while(!add.key(dec,bit,sum.one,F))
{
hash.level<-hash.level+1
bit1<- dectobit.alt(2654435761*(dec+97*sum.one+hash.level*36599)+hash.level*59+hash.level)
dec<- hashing(bit1)+1
#print(dec)
}
dec <- dec - 1
}
return(dec)
}
},
add.key = function(dec,bit,sum.one,levl)
{
lb<-length(bit)
if(dec>2^hash.length)
return(FALSE)
if(levl)
{
if(is.na(statistics1[dec,1]))
{
statistics1[dec,16]<-sum.one
hash.keys1[dec,]<-bit
return(TRUE)
}
if(is.na(statistics1[dec,16]))
{
statistics1[dec,16]<-sum.one
hash.keys1[dec,]<-bit#c(array(0,dim = Nvars-lb),bit)
return(TRUE)
}
if(statistics1[dec,16]!=sum.one)
return(FALSE)
i<-1
#print(lb)
lb <- length(bit)
while(i<=lb&&(hash.keys1[dec,Nvars-i+1])==bit[lb - i +1])
i<-i+1
return((i-1)==lb)
}
else
{
if(is.na(statistics[dec,1]))
{
statistics[dec,16]<-sum.one
hash.keys[dec,]<-bit
return(TRUE)
}
if(is.na(statistics[dec,16]))
{
statistics[dec,16]<-sum.one
hash.keys[dec,]<-bit#c(array(0,dim = Nvars-lb),bit)
return(TRUE)
}
if(statistics[dec,16]!=sum.one)
return(FALSE)
i<-1
lb <- length(bit)
while(i<=lb&&(hash.keys[dec,Nvars-i+1])==bit[lb - i +1])
i<-i+1
return((i-1)==lb)
}
},
hashing = function(bit)
{
n<-length(bit)
if(n<hash.length)
return(bittodec.alt(bit))
return(bittodec.alt(bit[(n-hash.length+1):n]))
},
binlog = function (x) # bitwise logorithm (base 2) computations based on Al Kashi's algorithm
{
lx <- length(x)
tol = -lx
y <- 0 # initialise output
b <- 0.5 # initialise mantissa
# arrange the input into a known range
if(lx == 1){
if(x[1] == 0)
return(-Inf)
if(x[1] == 1)
return(0)
}
# move one bit to the left end
# elsewhise
if(x[2]==0 && lx==2)
return(1)
powto<-2^(-c(1:(lx-1)))
float.x<-sum(x[2:lx]*powto)
x <- 1 + float.x
y <- lx -1
f<-0
fb<--1
# move one bit to the right end
# now x = 1.5
# loop until desired tolerance met
while(fb > tol)
{
x <- x*x
# update the index
if (x >= 2)
{
x <- x/2
y <- y + b
f <- log(exp(f)+exp(fb))
}
# scale for the next bit
b <- b/2
fb <- (fb-log(2))
}
return(list(y = y,z = lx -1, f = f))
},
#transform decimal numbers to binary
dectobit = function(dec) #transform a natural number into a binary vector to correspond between vector of solutions and storage array
{
if(!double.hashing || dec < 2^hash.length){
if(dec == 0)
return(0)
q<-dec
bin<-NULL
while(q!=0)
{
r<-q/2
q=floor(r)
bin<-c(as.integer(2*(r-q)),bin)
}
return(bin)
}
return(dehash(dec))
},
dectobit.alt = function(dec) #transform a natural number into a binary vector to correspond between vector of solutions and storage array
{
if(dec == 0)
return(0)
q<-dec
bin<-NULL
while(q!=0)
{
r<-q/2
q=floor(r)
bin<-c(as.integer(2*(r-q)),bin)
}
return(bin)
},
dehash=function(dec)
{
if(exists("hash.keys1"))
return(hash.keys1[dec,])
if(exists("hash.keys"))
return(hash.keys[dec,])
return(NULL)
},
#calculate move probabilities
calculate.move.logprobabilities = function(varold, varnew, switch.type,min.N,max.N)
{
if(switch.type == 1) # random size random N(x)
{
warning("This option should not be chosen for randomization unless p.add == 0.5 ", call. = FALSE)
min.N = max.N
log.mod.switch.prob <- log(1/(max.N - min.N +1)) # probability of having that many differences
KK<-sum(abs(varold-varnew))
log.mod.switch.prob <- 0 #always the same probabilities for moves within thenighbourhood for p=0.5
log.mod.switchback.prob <- 0
} else if(switch.type == 2) #fixed N(x) inverse operator
{
if(min.N!=max.N)
{
warning("min.N should be equal to max.N in swap type neighbourhoods min.N:=max.N", call. = FALSE)
min.N = max.N
}
log.mod.switch.prob <- log(1/(max.N - min.N +1))
KK<-max.N
log.mod.switch.prob <- log.mod.switch.prob + KK*log(factorial(Nvars - KK + 1)/factorial(Nvars))
log.mod.switchback.prob <- log.mod.switch.prob
}else if(switch.type == 3) # random sized inverse N(x)
{
log.mod.switch.prob <- log(1/(max.N - min.N +1))
KK<-sum(abs(varold-varnew))
log.mod.switch.prob <- log.mod.switch.prob + KK*log(factorial(Nvars - KK + 1)/factorial(Nvars))
log.mod.switchback.prob <- log.mod.switch.prob
}else if(switch.type == 4) # fixed N(x) for reverse from type 2 swaps
{
if(min.N!=max.N)
{
warning("min.N should be equal to max.N in swap type neighbourhoods min.N:=max.N", call. = FALSE)
min.N = max.N
}
log.mod.switch.prob <- log(1/(max.N - min.N +1))
KK<-max.N
log.mod.switch.prob <- log.mod.switch.prob + KK*log(factorial(Nvars - KK + 1)/factorial(Nvars))
log.mod.switchback.prob <- log.mod.switch.prob
}else if(switch.type > 4)
{
log.mod.switch.prob <- log(x = 1)
log.mod.switchback.prob <-log(x = 1)
}
return(list(log.switch.forw.prob = log.mod.switch.prob, log.switch.back.prob = log.mod.switchback.prob))
},
#build a new model to draw from the old one return selection probabilities
buildmodel= function(varcur.old,statid, shift.cpu,max.cpu,switch.type,min.N,max.N,changeble.coord = NULL)
{
vect<- vector(length = max.cpu,mode = "list")
shift<-0
for(cpu in 1:(max.cpu))
{
set.seed(runif(1,1,10000), kind = NULL, normal.kind = NULL)
if(!is.null(varcur.old))
{
varcur<-varcur.old
}else
{
varcur<-rbinom(n = (Nvars),size = 1,0.5)
varcur.old<-varcur
}
changeble <-FALSE
if(!is.null(changeble.coord))
{
changeble <- TRUE
}
if(switch.type == 1) # random size random N(x) # try avoiding when doing mcmc rather than optimization
{
log.mod.switch.prob <- log(1/(max.N - min.N +1))
KK<-floor(runif(n=1, min.N, max.N + 0.999999999))
log.mod.switch.prob <- log.mod.switch.prob + KK*log(factorial(Nvars - KK + 1)/factorial(Nvars))
log.mod.switchback.prob <- log.mod.switch.prob
change.buf <- array(data = 0,dim = Nvars)
if(changeble){
for(ttt in 1:KK)
{
iid<-floor(runif(n = 1,min = 1,max = Nvars+0.999999999))
if(changeble.coord[iid]==1)
next
change.buf[iid] = 1
varcur[iid] <-rbinom(n = 1,size = 1,prob = p.add[iid])
log.mod.switch.prob = 0#log.mod.switch.prob + log(dbinom(x = varcur[iid],size = 1,prob = p.add[iid])) # this is one of the pathes to get there in general
log.mod.switchback.prob = 0# log.mod.switchback.prob + log(dbinom(x = 1-varcur[iid],size = 1,prob = p.add[iid])) # but this is one of the ways to get back only
}
}else
{
iid<-floor(runif(n = KK,min = 1,max = Nvars+0.999999999))
change.buf[iid] = 1
varcur[iid] <-rbinom(n = KK,size = 1,prob = p.add)
}
}else if(switch.type == 2) # fixed sized inverse N(x)
{
if(min.N!=max.N)
{
warning("min.N should be equal to max.N in swap type neighbourhoods min.N:=max.N", call. = FALSE)
min.N <- max.N
}
log.mod.switch.prob <- log(1/(max.N - min.N +1))
KK<-max.N
log.mod.switch.prob <- log.mod.switch.prob + KK*log(factorial(Nvars - KK + 1)/factorial(Nvars))
log.mod.switchback.prob <- log.mod.switch.prob
change.buf <- array(data = 0,dim = Nvars)
if(changeble){
for(ttt in 1:KK)
{
iid<-floor(runif(n = 1,min = 1,max = Nvars+0.999999999))
if(change.buf[iid]==1 || changeble.coord[iid]==1)
{
KK<-KK+1
next
}
change.buf[iid] = 1
varcur[iid] = 1-varcur[iid]
}
}else
{
iid<-floor(runif(n = KK,min = 1,max = Nvars+0.999999999))
change.buf[iid] = 1
varcur[iid] = 1-varcur[iid]
}
}else if(switch.type == 3) # random sized inverse N(x)
{
log.mod.switch.prob <- log(1/(max.N - min.N +1))
KK<-floor(runif(n=1, min.N, max.N + 0.999999999))
log.mod.switch.prob <- log.mod.switch.prob + KK*log(factorial(Nvars - KK + 1)/factorial(Nvars))
log.mod.switchback.prob <- log.mod.switch.prob
change.buf <- array(data = 0,dim = Nvars)
if(changeble){
for(ttt in 1:KK)
{
iid<-floor(runif(n = 1,min = 1,max = Nvars+0.999999999))
if(change.buf[iid]==1 || changeble.coord[iid]==1)
{
KK<-KK+1
next
}
change.buf[iid] = 1
varcur[iid] = 1-varcur[iid]
}
}else
{
iid<-floor(runif(n = KK,min = 1,max = Nvars+0.999999999))
change.buf[iid] = 1
varcur[iid] = 1-varcur[iid]
}
}else if(switch.type == 4) # fixed N(x) for reverse from type 2 swaps
{
if(min.N!=max.N)
{
warning("min.N should be equal to max.N in swap type neighbourhoods min.N:=max.N", call. = FALSE)
min.N <- max.N
}
log.mod.switch.prob <- log(1/(max.N - min.N +1))
KK<-max.N
log.mod.switch.prob <- log.mod.switch.prob + KK*log(factorial(Nvars - KK + 1)/factorial(Nvars))
log.mod.switchback.prob <- log.mod.switch.prob
ids <- which(changeble.coord == 0)
change.buf <- array(data = 0,dim = Nvars)
if(changeble){
for(ttt in KK)
{
iid<-floor(runif(n = 1,min = 1,max = Nvars+0.999999999))
if(change.buf[iid]==1 || changeble.coord[iid]==1)
{
KK<-KK+1
next
}
change.buf[iid] = 1
varcur[iid] = 1-varcur[iid]
}}
else
{
iid<-floor(runif(n = KK,min = 1,max = Nvars+0.999999999))
change.buf[iid] = 1
varcur[iid] = 1-varcur[iid]
}
}else if(switch.type == 5)
{
change.buf <- array(data = 0,dim = (Nvars))
log.mod.switch.prob<-0
log.mod.switchback.prob <-0
add<-0
if(cpu+shift>Nvars)
{
shift<-1 - cpu
varcur.old<-rep(0,times = Nvars)
}
if(varcur.old[cpu+shift]!=1)
{
varcur<-varcur.old
varcur[cpu+shift]<-1
}else
{
if(cpu+shift<Nvars)
shift<-shift+1
while(varcur.old[cpu+shift]==1)
{
shift<-shift+1
if(cpu+shift>=Nvars)
{
shift <- (Nvars - cpu)
break
}
}
varcur<-varcur.old
varcur[cpu+shift]<-1
}
}else if(switch.type == 6)
{
change.buf <- array(data = 0,dim = (Nvars))
log.mod.switch.prob<-0
log.mod.switchback.prob <-0
add<-0
if(cpu+shift>Nvars)
{
shift<-1 - cpu
varcur.old<-rep(1,times = Nvars)
}
if(varcur.old[cpu+shift]!=0)
{
varcur<-varcur.old
varcur[cpu+shift]<-0
}else
{
if(cpu+shift<Nvars)
shift<-shift+1
while(varcur.old[cpu+shift]==0)
{
shift<-shift+1
if(cpu+shift>=Nvars)
{
shift <- (Nvars - cpu)
break
}
}
varcur<-varcur.old
varcur[cpu+shift]<-0
}
}else if(switch.type == 7)
{
change.buf <- array(data = 0,dim = (Nvars))
log.mod.switch.prob<-0
log.mod.switchback.prob <-0
vec<-dectobit(cpu + shift.cpu)
varcur<-c(array(0,dim = (Nvars -length(vec))),vec) #issues here
}else if(switch.type == 8)
{
log.mod.switch.prob <-0
log.mod.switchback.prob <-0
varcur<-varcur.old
change.buf <- array(data = 0,dim = Nvars)
#if(printable.opt)print("type 8 invoked")
#if(printable.opt)print(varcur)
}else
{
log.mod.switch.prob <- 0
log.mod.switchback.prob <-0
change.buf <- array(data = 1,dim = Nvars)
varcur <-rbinom(n = Nvars,size = 1,prob = p.add)
#if(printable.opt)print("type 8 invoked")
#if(printable.opt)print(varcur)
}
# for(g in 1:max(isobsbinary))
# {
# if(length(varcur[which(isobsbinary == g && varcur %in% c(1,3))])==length(which(isobsbinary == g)))
# varcur[which(isobsbinary == g && varcur %in% c(1,3))[1]]=0
# if(length(varcur[which(isobsbinary == g && varcur %in% c(2,3))])==length(which(isobsbinary == g)))
# varcur[which(isobsbinary == g && varcur %in% c(2,3))[1]]=0
# }
covobs <- if(fparam[1]=="Const")fparam[which(varcur[-1] == 1)+1]else fparam[which(varcur==1)]
#obsconst<-2*as.integer((varcur[1] %in% c(1,3)) || length(covobs) == 0 || (length(covobs)==length(which(isobsbinary != 0))) ) -1
obsconst<-ifelse(fparam[1]=="Const",2*as.integer((varcur[1])) -1,1)
id<-bittodec(varcur)
id<-id+1
if(ifelse(exists("statistics1"),is.na(statistics1[id,1]),ifelse(exists("statistics"),is.na(statistics[id,1,]),ifelse(exists("hashStat"),!has.key(hash = hashStat,key = paste(varcur,collapse = "")),TRUE))))#||TRUE)
{
formula <- NULL
capture.output({withRestarts(tryCatch(capture.output({formula <- as.formula(paste(paste(fobserved[1]), " ~ ",obsconst,ifelse(length(covobs)>0," + ",""), paste(covobs, collapse=" + "), latent.formula)) })), abort = function(){onerr<-TRUE;fm<-NULL})}) ## not considered currently in RJMCMC, is only valid for model selection
if(is.null(formula)){
formula <- as.formula(paste(paste(fobserved[1]), " ~ ",obsconst,ifelse(length(covobs)>0," + ",""), paste(covobs, collapse=" + ")))
}
}else
{
formula <- NULL
}
vect[[cpu]]<-list(formula = formula, varcur = varcur, statid = statid, changed = change.buf, log.mod.switch.prob = log.mod.switch.prob , log.mod.switchback.prob = log.mod.switchback.prob )
}
return(vect)
},
#fit the given model by means of the specified estimator
fitmodel = function(model)
{
if(!is.null(model))
{
fm<-NULL
id<-bittodec(model$varcur)
if(is.null(id))
id = 0
id<-id+1
if(exists("statistics1")){
if(is.na(statistics1[id,1]))
{
#if(printable.opt)print("Invoked from EMJMCMC environment")
onerr<-FALSE
#if(printable.opt)print("INLA internal error")
statistics1[id,c(2,3)]<-100000
statistics1[id,1]<--100000
statistics1[id,4:14]<-0
#prior = "normal",param = c(model$beta.mu.prior,model$beta.tau.prior))
# capture.output({withRestarts(tryCatch(capture.output({fm<-inla(formula = model$formula,family = "binomial",Ntrials = data$total_bases,data = data,control.fixed = list(mean = list(default = model$beta.mu.prior),mean.intercept = model$beta.mu.prior, prec = list(default = model$beta.tau.prior), prec.intercept = model$beta.tau.prior) ,control.compute = list(dic = model$dic.t, waic = model$waic.t, mlik = model$mlik.t))
# })), abort = function(){onerr<-TRUE;fm<-NULL})}) # fit the modal, get local improvements
#
capture.output({withRestarts(tryCatch(capture.output({fm<-do.call(estimator, c(estimator.args, model$formula))
})), abort = function(){onerr<-TRUE;fm<-NULL})}) # fit the modal, get local improvements
if(!is.null(fm))
{
statistics1[id,2]<-fm$waic[[1]]
statistics1[id,1]<-fm$mlik[[1]]
statistics1[id,3]<-fm$dic[[1]]
# if(id>1)
# {
# inxx<-which(model$varcur==1)
# if(length(inxx)==length(fm$summary.fixed$mean))
# statistics1[id,15+inxx]<-fm$summary.fixed$mean
# }
if(fm$waic[[1]]<g.results[2,1] && !is.na(fm$waic[[1]]))
{
g.results[2,1]<-fm$waic[[1]]
g.results[2,2]<-as.integer(id)
}
if(fm$mlik[[1]]>g.results[1,1] && !is.na(fm$mlik[[1]]))
{
g.results[1,1]<-fm$mlik[[1]]
g.results[1,2]<-as.integer(id)
}
if(fm$dic[[1]]<g.results[3,1]&& !is.na(fm$dic[[1]]))
{
g.results[3,1]<-fm$dic[[1]]
g.results[3,2]<-as.integer(id)
}
g.results[4,2] <- g.results[4,2]+1
if(g.results[4,2]%%recalc.margin == 0)
{
p.add <<- as.array(post_proceed_results(statistics1)$p.post)
}
}
}
if(model$statid!=-1)
statistics1[id,model$statid+1]<-statistics1[id,model$statid+1] + 1
g.results[4,1] <- g.results[4,1]+1
return(list(mlik=statistics1[id,1],waic=statistics1[id,2],dic=statistics1[id,3]))
}else if(exists("statistics")){
if(is.na(statistics[id,1]))
{
#if(printable.opt)print("Invoked from EMJMCMC SUB environment")
onerr<-FALSE
#if(printable.opt)print("INLA internal error")
statistics[id,c(2,3)]<-100000
statistics[id,1]<--100000
statistics[id,4:14]<-0
#prior = "normal",param = c(model$beta.mu.prior,model$beta.tau.prior))
# capture.output({withRestarts(tryCatch(capture.output({fm<-inla(formula = model$formula,family = "binomial",Ntrials = data$total_bases,data = data,control.fixed = list(mean = list(default = model$beta.mu.prior),mean.intercept = model$beta.mu.prior, prec = list(default = model$beta.tau.prior), prec.intercept = model$beta.tau.prior) ,control.compute = list(dic = model$dic.t, waic = model$waic.t, mlik = model$mlik.t))
# })), abort = function(){onerr<-TRUE;fm<-NULL})}) # fit the modal, get local improvements
#
capture.output({withRestarts(tryCatch(capture.output({fm<-do.call(estimator, c(estimator.args, model$formula))
})), abort = function(){onerr<-TRUE;fm<-NULL})}) # fit the modal, get local improvements
if(!is.null(fm))
{
statistics[id,2]<-fm$waic[[1]]
statistics[id,1]<-fm$mlik[[1]]
statistics[id,3]<-fm$dic[[1]]
# if(id>1)
# {
# inxx<-which(model$varcur==1)
# if(length(inxx)==length(fm$summary.fixed$mean))
# statistics[id,14+inxx]<-fm$summary.fixed$mean
# }
if(fm$waic[[1]]<g.results[2,1] && !is.na(fm$waic[[1]]))
{
g.results[2,1]<-fm$waic[[1]]
g.results[2,2]<-as.integer(id)
}
if(fm$mlik[[1]]>g.results[1,1] && !is.na(fm$mlik[[1]]))
{
g.results[1,1]<-fm$mlik[[1]]
g.results[1,2]<-as.integer(id)
}
if(fm$dic[[1]]<g.results[3,1]&& !is.na(fm$dic[[1]]))
{
g.results[3,1]<-fm$dic[[1]]
g.results[3,2]<-as.integer(id)
}
g.results[4,2] <- g.results[4,2]+1
if(g.results[4,2]%%recalc.margin == 0)
{
proceeeded <- post_proceed_results(statistics)
p.add <<- as.array(proceeeded$p.post)
#g.results[4,2] <-
}
}
}
if(model$statid!=-1)
statistics[id,model$statid+1]<-statistics[id,model$statid+1] + 1
g.results[4,1] <- g.results[4,1]+1
return(list(mlik=statistics[id,1],waic=statistics[id,2],dic=statistics[id,3]))
}else if(exists("hashStat")){
idd<- as.character(paste(model$varcur,collapse = ""))
if(!has.key(key = idd,hash = hashStat))
{
#if(printable.opt)print("Invoked from EMJMCMC hash table environment")
onerr<-FALSE
#if(printable.opt)print("INLA internal error")
hashBuf<-array(data = 0,dim = 3)
hashBuf[c(2,3)]<-100000
hashBuf[1]<- -100000
#prior = "normal",param = c(model$beta.mu.prior,model$beta.tau.prior))
# capture.output({withRestarts(tryCatch(capture.output({fm<-inla(formula = model$formula,family = "binomial",Ntrials = data$total_bases,data = data,control.fixed = list(mean = list(default = model$beta.mu.prior),mean.intercept = model$beta.mu.prior, prec = list(default = model$beta.tau.prior), prec.intercept = model$beta.tau.prior) ,control.compute = list(dic = model$dic.t, waic = model$waic.t, mlik = model$mlik.t))
# })), abort = function(){onerr<-TRUE;fm<-NULL})}) # fit the modal, get local improvements
#
capture.output({withRestarts(tryCatch(capture.output({fm<-do.call(estimator, c(estimator.args, model$formula))
})), abort = function(){onerr<-TRUE;fm<-NULL})}) # fit the modal, get local improvements
if(!is.null(fm))
{
hashBuf[2]<-fm$waic[[1]]
hashBuf[1]<-fm$mlik[[1]]
hashBuf[3]<-fm$dic[[1]]
hashStat[idd] <- hashBuf
# if(id>1)
# {
# inxx<-which(model$varcur==1)
# if(length(inxx)==length(fm$summary.fixed$mean))
# statistics[id,14+inxx]<-fm$summary.fixed$mean
# }
if(fm$waic[[1]]<g.results[2,1] && !is.na(fm$waic[[1]]))
{
g.results[2,1]<-fm$waic[[1]]
g.results[2,2]<-(id)
}
if(fm$mlik[[1]]>g.results[1,1] && !is.na(fm$mlik[[1]]))
{
g.results[1,1]<-fm$mlik[[1]]
g.results[1,2]<-(id)
}
if(fm$dic[[1]]<g.results[3,1]&& !is.na(fm$dic[[1]]))
{
g.results[3,1]<-fm$dic[[1]]
g.results[3,2]<-(id)
}
g.results[4,2] <- g.results[4,2]+1
if(g.results[4,2]%%recalc.margin == 0)
{
proceeeded <- post_proceed_results_hash(hashStat)
p.add <<- as.array(proceeeded$p.post)
#g.results[4,2] <-
}
}
}
g.results[4,1] <- g.results[4,1]+1
hasRes<- values(hashStat[idd])
return(list(mlik=hasRes[1],waic=hasRes[2],dic=hasRes[3]))
}else
{
capture.output({withRestarts(tryCatch(capture.output({fm<-do.call(estimator, c(estimator.args, model$formula))
})), abort = function(){onerr<-TRUE;fm<-NULL})}) # fit the modal, get local improvements
if(!is.null(fm)){
if(fm$waic[[1]]<g.results[2,1] && !is.na(fm$waic[[1]]))
{
g.results[2,1]<-fm$waic[[1]]
g.results[2,2]<-(id)
}
if(fm$mlik[[1]]>g.results[1,1] && !is.na(fm$mlik[[1]]))
{
g.results[1,1]<-fm$mlik[[1]]
g.results[1,2]<-(id)
}
if(fm$dic[[1]]<g.results[3,1]&& !is.na(fm$dic[[1]]))
{
g.results[3,1]<-fm$dic[[1]]
g.results[3,2]<-(id)
}
g.results[4,1] <- g.results[4,1]+1
g.results[4,2] <- g.results[4,2]+1
return(list(mlik=fm$mlik[[1]],waic=fm$waic[[1]],dic=fm$dic[[1]]))
}
else
{
g.results[4,1] <- g.results[4,1]+1
return(list(mlik=-Inf,waic=Inf,dic=Inf))
}
}
}
g.results[4,1] <- g.results[4,1]+1
g.results[4,2] <- g.results[4,2]+1
return(list(mlik=-Inf,waic=Inf,dic=Inf))
},
#lambda function for mtmcmc
lambda = function(c,alpha,g1,g2,g.domain.pos) # simmetric choice driving function
{
if((c!=0))
{
res<-((((1/c)*(1+(g1+g2)*(g.domain.pos))/(1-(g1+g2)*(1-g.domain.pos)))^alpha))
#!#if(printable.opt)print(res)
return(res)
}else
{
#!#if(printable.opt)print(1)
return(1)
}
},
#norm between probabilities vectors
normprob=function(p1,p2)
{
nn<-abs(1-sum((p1+0.1)/(p2+0.1))/length(p1))
if(is.na(nn))
nn<-Inf
return(nn)
},
#local mcmc procedure
learnlocalMCMC=function(model)
{
M<-M.mcmc
mlikcur<- model$mlikcur
waiccur<-model$waiccur
varcand<-model$varcur
varcur<-model$varcur
varcurb<-model$varcur
varglob<-model$varcur
modglob<-NULL
fm<-NULL
fmb<-NULL
# estimate large jump in a reverse move
if(model$reverse || is.infinite(mlikcur))
{
vectbg<-buildmodel(max.cpu = 1,varcur.old = varcur,statid = model$statid,switch.type=8, min.N = min.N,max.N = max.N)
if(!is.null(vectbg[[1]]$formula))
{
bgmod <- lapply(X = vectbg,FUN = .self$fitmodel)
waiccur<-bgmod[[1]]$waic
mlikcur<-bgmod[[1]]$mlik
}
else if(exists("statistics1"))
{
iidd<-bittodec(varcur)+1
waiccur<-statistics1[iidd,2]
mlikcur<-statistics1[iidd,1]
}else if(exists("hashStat"))
{
iidd<-paste(varcur,collapse = "")
waiccur<-values(hashStat[iidd])[2]
mlikcur<-values(hashStat[iidd])[1]
}
}
if(printable.opt)print(paste("Begin with ", mlikcur))
mlikglob<- mlikcur
mlikcand<- mlikcur
waiccand<- waiccur
waicglob<- waiccur
waiccur<- waiccur
for(m in 1:M)
{
withRestarts(tryCatch({
# statistics <- describe(statistics)
vect<-buildmodel(varcur.old = varcur,statid = model$statid,max.cpu = max.cpu,switch.type=switch.type, min.N = min.N, max.N = max.N)
cluster<-TRUE
flag1<-0
for(mod_id in 1:max.cpu)
{
if(is.null(vect[[mod_id]]$formula))
{
flag1<-flag1+1
}
}
#flag1<-sum(is.null(vect[[]]$formula))
if(flag1==max.cpu)
{
cluster<-FALSE
if(printable.opt)print("!!!!MTMCMC models already estimated!!!!")
}else
{
res.par <- parallelize(X = vect,FUN = .self$fitmodel)
}
p.select.y <- array(data = 0, dim = max.cpu)
for(mod_id in 1:max.cpu)
{
if(cluster)
{
fm<-res.par[[mod_id]]
if(is.null(fm)&&(is.na(res.par[[mod_id]]$waic)))
{
varcand<-varcurb
if(printable.opt)print("locMTMCMC Model Fit Error!?")
next
}
}
varcand<-vect[[mod_id]]$varcur
if(cluster)
{
waiccand<-res.par[[mod_id]]$waic
mlikcand<-res.par[[mod_id]]$mlik
}else if(exists("statistics1"))
{
iidd<-bittodec(varcand)+1
waiccand<-statistics1[iidd,2]
mlikcand<-statistics1[iidd,1]
}else if(exists("hashStat"))
{
iidd<-paste(varcand,collapse = "")
waiccand<-values(hashStat[iidd])[2]
mlikcand<-values(hashStat[iidd])[1]
}
if((mlikcand>mlikglob)) #update the parameter of interest
{
if(printable.opt)print(paste("locMTMCMC update waic.glob = ", waiccand))
if(printable.opt)print(paste("locMTMCMC update waic.glob.mlik = ", mlikglob))
mlikglob<-mlikcand
waicglob<-waiccand
varglob<-varcand
if(cluster)
modglob<-fm
}
g1 <- waiccur
if(waiccur == Inf)
{
g1 = 1
}
p.select.y[mod_id]<-(mlikcand + vect[[mod_id]]$log.mod.switchback.prob+log(lambda(c = cc, alpha = aa, g1 = -g1, g2 = -waiccand,g.domain.pos = FALSE))) # correct for different criteria later
if(is.na(p.select.y[mod_id]))
p.select.y[mod_id] <- 0
if(is.infinite(p.select.y[mod_id]) || p.select.y[mod_id]>100000000)
{
#if(printable.opt)print(paste("very large log.w.y detected ",p.select.y[mod_id]))
p.select.y[mod_id] <- 100000000
}
}
max.p.select.y <- max(p.select.y)
p.select.y<-p.select.y-max.p.select.y
#if(printable.opt)print(paste("max log.w.y is ",max.p.select.y,"normilized log.w.n.y is ", paste(p.select.y,collapse = ", ")))
ID<-sample(x = max.cpu,size = 1,prob = exp(p.select.y))
if(printable.opt)print(paste("cand ",ID," selected"))
varcand<-vect[[ID]]$varcur
if(cluster)
{
waiccand<-res.par[[ID]]$waic
mlikcand<-res.par[[ID]]$mlik
}else if(exists("statistics1"))
{
iidd<-bittodec(varcand)+1
waiccand<-statistics1[iidd,2]
mlikcand<-statistics1[iidd,1]
}else if(exists("hashStat"))
{
iidd<-paste(varcand,collapse = "")
waiccand<-values(hashStat[iidd])[2]
mlikcand<-values(hashStat[iidd])[1]
}
#p.Q.cand<- p.select.y[ID]/sum(p.select.y)
if(printable.opt)print("do reverse step")
p.select.z <- array(data = 0.01, dim = max.cpu)
if(max.cpu!=1)
{
vect1<-buildmodel(max.cpu = max.cpu -1,varcur.old = varcand,statid = model$statid,switch.type = switch.type, min.N = min.N, max.N = max.N)
cluster<-TRUE
flag1<-0
for(mod_id in 1:(max.cpu-1))
{
if(is.null(vect1[[mod_id]]$formula))
{
flag1<-flag1+1
}
}
if(flag1==(max.cpu-1))
{
cluster<-FALSE
if(printable.opt)print("!!!!MTMCMC reverse models already estimated!!!!")
}else
{
res.par.back <- parallelize(X = vect1,FUN = .self$fitmodel)
}
for(mod_id in 1:(max.cpu-1))
{
if(cluster)
{
if(is.null(fm)&&(is.na(res.par.back[[mod_id]]$waic)))
{
if(printable.opt)print("locMTMCMC Model Fit Error!?")
next
}
}
varcand.b<-vect1[[mod_id]]$varcur
if(cluster)
{
waiccand.b<-res.par.back[[mod_id]]$waic
mlikcand.b<-res.par.back[[mod_id]]$mlik
}else if(exists("statistics1"))
{
iidd<-bittodec(varcand.b)+1
waiccand.b<-statistics1[iidd,2]
mlikcand.b<-statistics1[iidd,1]
}else if(exists("hashStat"))
{
iidd<-paste(varcand.b,collapse = "")
waiccand.b<-values(hashStat[iidd])[2]
mlikcand.b<-values(hashStat[iidd])[1]
}
if((mlikcand.b>mlikglob))
{
if(printable.opt)print(paste("locMTMCMC update waic.glob = ", waiccand.b))
if(printable.opt)print(paste("locMTMCMC update waic.glob.mlik = ", mlikcand.b))
mlikglob<-mlikcand.b
waicglob<-waiccand.b
varglob<-varcand.b
if(cluster)
modglob<-fm
}
g1 = waiccand
if(waiccand == Inf)
{
g1 = 1
}
p.select.z[mod_id]<-(mlikcand.b+vect1[[mod_id]]$log.mod.switchback.prob+(lambda(c = cc, alpha = aa, g1 = -g1, g2 = -waiccand.b,g.domain.pos = FALSE))) # correct for different criteria later
if(is.na(p.select.z[mod_id]))
p.select.z[mod_id]=0
if(is.infinite(p.select.z[mod_id]) || p.select.z[mod_id] > 100000000)
{
#if(printable.opt)print(paste("very large log.w.y detected ",p.select.z[mod_id]))
p.select.z[mod_id] <- 100000000
}
}
}
if( waiccur == Inf)
{
g1 = 1
}
p.select.z[max.cpu] <- (mlikcur+vect[[ID]]$log.mod.switch.prob+(lambda(c = cc, alpha = aa, g1 = -g1, g2 = -waiccand,g.domain.pos = FALSE)))
if(is.na(p.select.z[mod_id]))
p.select.z[mod_id]=0
if(is.infinite(p.select.z[mod_id]) || p.select.z[mod_id] > 100000000)
{
#if(printable.opt)print(paste("very large log.w.y detected ",p.select.z[mod_id]))
p.select.z[mod_id] <- 100000000
}
max.p.select.z <- max(p.select.z)
p.select.z<-p.select.z-max.p.select.z
if(printable.opt)print(paste("max log.w.z is ",max.p.select.z,"normilized log.w.n.z is ", paste(p.select.z,collapse = ", ")))
if(log(runif(n = 1,min = 0,max = 1)) < (log(sum(exp(p.select.y)))-log(sum(exp(p.select.z)))) + max.p.select.y - max.p.select.z )
{
mlikcur<-mlikcand
if(printable.opt)print(paste("locMTMCMC update ratcur = ", mlikcand))
if(printable.opt)print(paste("locMTMCMC accept move with ", waiccand))
varcur<-varcand
waiccur<-waiccand
}
}),abort = function(){fm<-fmb;closeAllConnections();options(error=traceback); onerr<-TRUE})
}
#if(printable.opt)print("FINISH LOCAL MTMCMC")
#!#if(printable.opt)print(points)
if(model$reverse == FALSE)
{
if(is.null(varcur))
{
#if(printable.opt)print("No moves acceoted in the procedure")
varcur<-varcand
waiccur<-waiccand
varcur<-varcand
}
vect<-buildmodel(max.cpu = 1,varcur.old = varcur,statid = model$statid,switch.type = type.randomize,min.N = min.N.randomize,max.N = max.N.randomize)
varcur<-vect[[1]]$varcur
#if(printable.opt)print(varcur)
cluster<-TRUE
if(is.null(vect[[1]]$formula))
{
cluster<-FALSE
if(printable.opt)print("!!!!MTMCMC reverse models already estimated!!!!")
}else
{
mod<-fitmodel(vect[[1]])
}
if(cluster)
{
waiccur<-mod$waic
mlikcur<-mod$mlik
}else if(exists("statistics1"))
{
iidd<-bittodec(varcur)+1
waiccur<-statistics1[iidd,2]
mlikcur<-statistics1[iidd,1]
}else if(exists("hashStat"))
{
iidd<-paste(varcur,collapse = "")
waiccur<-values(hashStat[iidd])[2]
mlikcur<-values(hashStat[iidd])[1]
}
# incorporate what happens for the backward optimization
model.prob<-vect[[1]]$log.mod.switch.prob
model.prob.fix<-vect[[1]]$log.mod.switchback.prob
}else # incorporate what happens for the reverse move
{
if(is.null(varcur))
{
#if(printable.opt)print("No moves acceoted in the reverse procedure")
varcur<-varcand
waiccur<-waiccand
}
model.probs<-calculate.move.logprobabilities(varold = varcur, varnew = model$varold,switch.type = type.randomize,min.N = min.N.randomize,max.N = max.N.randomize)
model.prob<-model.probs$log.switch.forw.prob
model.prob.fix<-model.probs$log.switch.back.prob
}
if(is.null(varcur))
{
#if(printable.opt)print("NO VARCUR OBTAINED")
varcur<-i
waiccur<-Inf
varcur<-model$varold
}
return(list(varcur = varcur, waiccur = waiccur, mlikcur = mlikcur, log.prob.cur = model.prob,log.prob.fix = model.prob.fix, varglob = varglob, waicglob = waicglob, mlikglob = mlikglob, modglob = modglob))
},
#local simulated annealing optimization
learnlocalSA=function(model)
{
t.min<-SA.param$t.min
t<-SA.param$t.init
dt<-SA.param$dt
M<-SA.param$M
varcand<-model$varcur
varcurb<-model$varcur
varcur<-model$varcur
varglob<-model$varcur
varglob<-NULL
modglob<-NULL
probcur<-1
probrevcur<-1
first.prob <- 1
fm<-NULL
fmb<-NULL
mlikcur<- model$mlikcur
waiccur<-model$waiccur
# estimate large jump in a reverse move
# estimate large jump in a reverse move
if((model$reverse && !model$sa2) || is.infinite(mlikcur))
{
vectbg<-buildmodel(max.cpu = 1,varcur.old = varcur,statid = model$statid,switch.type=8, min.N = min.N,max.N = max.N)
if(!is.null(vectbg[[1]]$formula))
{
bgmod <- lapply(X = vectbg,FUN = .self$fitmodel)
waiccur<-bgmod[[1]]$waic
mlikcur<-bgmod[[1]]$mlik
}
else if(exists("statistics1"))
{
iidd<-bittodec(varcur)+1
waiccur<-statistics1[iidd,2]
mlikcur<-statistics1[iidd,1]
}else if(exists("hashStat"))
{
iidd<-paste(varcur,collapse = "")
waiccur<-values(hashStat[iidd])[2]
mlikcur<-values(hashStat[iidd])[1]
}
}
if(printable.opt)print(paste("Begin with ", mlikcur))
mlikglob<- mlikcur
mlikcand<- mlikcur
waiccand<- waiccur
waicglob<- waiccur
waiccur<- waiccur
while(t>t.min)
{
if(printable.opt)print(paste("anneal to ",t))
t.new<-t*exp(-dt)
if(model$reverse == TRUE && t.new <= t.min)
{
M<-M -1
}
for(m in 1:M)
{
withRestarts(tryCatch({
mmax.cpu = max.cpu
if(model$switch.type == 5)
{
mmax.cpu = Nvars - sum(varcur)
}else if(model$switch.type == 6)
{
mmax.cpu = sum(varcur)
}
if(mmax.cpu == 0)
mmax.cpu = 1
vect<-buildmodel(max.cpu = mmax.cpu,varcur.old = varcur,statid = model$statid,switch.type=model$switch.type, min.N = min.N,max.N = max.N)
cluster<-TRUE
flag1<-0
for(mod_id in 1:mmax.cpu)
{
if(is.null(vect[[mod_id]]$formula))
{
flag1<-flag1+1
}
}
if(flag1==mmax.cpu)
{
cluster<-FALSE
if(printable.opt)print("!!!!SA Models already estimated!!!!")
}else
{
res.par <- parallelize(X = vect,FUN = .self$fitmodel)
}
for(mod_id in 1:mmax.cpu)
{
if(cluster)
{
fmb<-fm
fm<-res.par[[mod_id]]
waiccand<-Inf
if(is.null(fm)&&(is.na(res.par[[mod_id]]$waic)))
{
varcand<-varcurb
if(printable.opt)print("SA Model Fit Error!?")
next
}
}
varcand<-vect[[mod_id]]$varcur
if(cluster)
{
waiccand<-res.par[[mod_id]]$waic
mlikcand<-res.par[[mod_id]]$mlik
}else if(exists("statistics1"))
{
iidd<-bittodec(varcand)+1
waiccand<-statistics1[iidd,2]
mlikcand<-statistics1[iidd,1]
}else if(exists("hashStat"))
{
iidd<-paste(varcand,collapse = "")
waiccand<-values(hashStat[iidd])[2]
mlikcand<-values(hashStat[iidd])[1]
}
if(objective==0)
{
objcand<-waiccand
objcur<-waiccur
objglob<-waicglob
}else
{
objcand<- -mlikcand
objcur<- -mlikcur
objglob<- -mlikglob
}
if(t == SA.param$t.init && mod_id == 1 && m == 2)
{
delta<-objcand - objcur
first.prob <- vect[[mod_id]]$log.mod.switchback.prob + log(punif(q = exp(x = delta/t),min = 0,max = 1))
}
if(objcand<objcur)
{
if(printable.opt)print(paste("SA accept move with ", objcand))
waiccur<-waiccand
varcur<-varcand
mlikcur<-mlikcand
if(objcand<objglob)
{
waicglob<-waiccand
varglob<-varcand
mlikglob<-mlikcand
if(cluster)
modglob<-fm
if(printable.opt)print(paste("SA update global optima with", objcand))
}
}else
{
delta<-objcand - objcur
if(runif(n = 1,min = 0,max = 1) <= exp(x = -delta/t))
{
model.probs<-calculate.move.logprobabilities(varold = varcur, varnew = varcand,switch.type = model$switch.type,min.N = min.N,max.N = max.N)
probcur<-model.probs$log.switch.forw.prob
probrevcur<-model.probs$log.switch.back.prob
waiccur<-waiccand
varcur<-varcand
mlikcur<-mlikcand
if(printable.opt)print(paste("SA accept move with ", objcand))
}
}
}
}),abort = function(){varcur<-varcurb; fm<-fmb;closeAllConnections();options(error=traceback); onerr<-TRUE})
}
t<-t.new
}
t<-t/exp(-dt)
if(model$reverse == FALSE)
{
model.prob<-log(punif(q = exp(x = -delta/t),min = 0,max = 1)) + probcur # log(P(Mk,Mk-1))
model.prob.fix<-log(punif(q = exp(x = delta/t),min = 0,max = 1)) + probrevcur # log(P(Mk-1,Mk))
if(model$sa2 == TRUE)
{
model.prob.fix<-model.prob.fix + first.prob # correcting for the term for local improvements of type 3.
}
}else # incorporate what happens for the reverse move
{
if(is.null(varcur))
{
if(printable.opt)print("No moves accepted in the reverse procedure")
varcur<-model$varcur
objcur<-model$objcur
}
delta<-objcur-model$objold
model.probs<-calculate.move.logprobabilities(varold = varcur, varnew = model$varold, switch.type = model$switch.type,min.N = min.N,max.N = max.N)
model.prob<-punif(q = exp(x = -delta/t),min = 0,max = 1,log.p = TRUE) + model.probs$log.switch.forw.prob
model.prob.fix<-punif(q = exp(x = delta/t),min = 0,max = 1,log.p = TRUE) + model.probs$log.switch.back.prob
if(model.prob==-Inf)
{
model.prob<- -100000000
}
if(model.prob.fix==-Inf)
{
model.prob.fix<- -100000000
}
}
return(list(varcur = varcur, waiccur = waiccur, mlikcur = mlikcur, log.prob.cur = model.prob,log.prob.fix = model.prob.fix, varglob = varglob, waicglob = waicglob, mlikglob = mlikglob, modglob = modglob))
},
#forward selection procedure
forward_selection=function(model)
{
if(printable.opt)print("begin forward selection procedure")
varcand<-model$varcur
varcurb<-model$varcur
varglob<-varcand
mlikglob<- model$mlikcur
mlikcur<- model$mlikcur
waiccand<- model$waiccur
waicglob<-model$waiccur
waiccur<-model$waiccur
waiccurb<-model$waiccur
varglob<-NULL
modglob<-NULL
fm<-NULL
fmb<-NULL
ub<-bittodec(array(1,length(varcurb)))
layer<-length(which(varcurb == 0))
mlikcur<- model$mlikcur
waiccur<-model$waiccur
# estimate large jump in a reverse move
# estimate large jump in a reverse move
if(is.infinite(mlikcur))
{
vectbg<-buildmodel(max.cpu = 1,varcur.old = varcurb,statid = model$statid,switch.type=8, min.N = min.N,max.N = max.N)
if(!is.null(vectbg[[1]]$formula))
{
bgmod <- lapply(X = vectbg,FUN = .self$fitmodel)
waiccur<-bgmod[[1]]$waic
mlikcur<-bgmod[[1]]$mlik
}
else if(exists("statistics1"))
{
iidd<-bittodec(varcand)+1
waiccur<-statistics1[iidd,2]
mlikcur<-statistics1[iidd,1]
}
else if(exists("hashStat"))
{
iidd<-paste(varcand,collapse = "")
waiccur<-values(hashStat[iidd])[2]
mlikcur<-values(hashStat[iidd])[1]
}
if(!is.na(mlikcur) && !is.na(waiccur) )
{
mlikglob<- mlikcur
mlikcur<- mlikcur
waiccand<- waiccur
waicglob<- waiccur
waiccur<- waiccur
waiccurb<- waiccur
}
}
while(layer>0)
{
withRestarts(tryCatch({
if(printable.opt)print(paste("proceed with layer",layer))
if(printable.opt)print(paste("current solution is",as.character(varcand)))
vect<-buildmodel(max.cpu = layer,varcur.old = varcurb,statid = model$statid, switch.type = 5,min.N = min.N, max.N = max.N)
if(printable.opt)print(paste("finish preparing models at layer",layer))
cluster<-TRUE
flag1<-0
for(mod_id in 1:layer)
{
if(is.null(vect[[mod_id]]$formula))
{
flag1<-flag1+1
}
}
if(flag1==layer)
{
cluster<-FALSE
if(printable.opt)print("!!!!forward Models already estimated!!!!")
}else
{
res.par <- parallelize(X = vect,FUN = .self$fitmodel)
}
if(printable.opt)print(paste("end forward optimizing at layer",layer))
for(mod_id in 1:layer)
{
if(cluster){
fmb<-fm
fm<-res.par[[mod_id]]
waiccand<-Inf
if(is.null(fm)&&(is.na(res.par[[mod_id]]$waic)))
{
varcand<-varcurb
if(printable.opt)print("forward Model Fit Error!?")
next
}
}
varcand<-vect[[mod_id]]$varcur
if(cluster)
{
waiccand<-res.par[[mod_id]]$waic
mlikcand<-res.par[[mod_id]]$mlik
}else if(exists("statistics1"))
{
iidd<-bittodec(varcand)+1
waiccand<-statistics1[iidd,2]
mlikcand<-statistics1[iidd,1]
}else if(exists("hashStat"))
{
iidd<-paste(varcand,collapse = "")
waiccand<-values(hashStat[iidd])[2]
mlikcand<-values(hashStat[iidd])[1]
}
if(objective==0)
{
objcand<-waiccand
objcur<-waiccur
objglob<-waicglob
}else
{
objcand<- -mlikcand
objcur<- -mlikcur
objglob<- -mlikglob
}
if(objcand<=objcur || mod_id ==1)
{
if(printable.opt)print(paste("forward accept with ", objcand))
objcur<-objcand
varcurb<-varcand
waiccur<-waiccand
varcur<-varcand
mlikcur<-mlikcand
if(objcur<objglob)
{
objglob<-objcur
waicglob<-waiccand
varglob<-varcand
mlikglob<-mlikcand
if(!is.null(fm))
modglob<-fm
if(printable.opt)print(paste("forward global optima with ", objcand))
}
}
# else
# {
# if(waiccand<=waiccur)
# {
# waiccur<-waiccand
# varcur<-varcand
# }
#
# }
}
if(objcur!=objglob)
{
if(model$locstop)
{
break
}else
{
varcurb<-varcur
}
}
}),abort = function(){varcur<-varcurb; fm<-fmb;closeAllConnections();options(error=traceback); onerr<-TRUE})
layer<-layer-1
}
model.prob<-1
model.prob.fix<-1
return(list(varcur = varglob, waiccur = waicglob, mlikcur = mlikglob, log.prob.cur = model.prob,log.prob.fix = model.prob.fix, varglob = varglob, waicglob = waicglob, mlikglob = mlikglob, modglob = modglob))
},
#backward selection procedure
backward_selection=function(model)
{
#if(printable.opt)print("begin backward selection procedure")
varcand<-model$varcur
varcurb<-model$varcur
varglob<-varcand
mlikglob<- model$mlikcur
mlikcur<- model$mlikcur
waiccur<- model$waiccur
waicglob<-model$waiccur
varglob<-NULL
modglob<-NULL
waiccurb<- model$waiccur
fm<-NULL
fmb<-NULL
ub<-bittodec(array(1,length(varcurb)))
layer<-length(which(varcurb == 1))
if(is.infinite(mlikcur))
{
vectbg<-buildmodel(max.cpu = 1,varcur.old = varcurb,statid = model$statid,switch.type=8, min.N = min.N,max.N = max.N)
if(!is.null(vectbg[[1]]$formula))
{
bgmod <- lapply(X = vectbg,FUN = .self$fitmodel)
waiccur<-bgmod[[1]]$waic
mlikcur<-bgmod[[1]]$mlik
}
else if(exists("statistics1"))
{
iidd<-bittodec(varcand)+1
waiccur<-statistics1[iidd,2]
mlikcur<-statistics1[iidd,1]
}else if(exists("hashStat"))
{
iidd<-paste(varcand,collapse = "")
waiccur<-values(hashStat[iidd])[2]
mlikcur<-values(hashStat[iidd])[1]
}
if(!is.na(mlikcur) && !is.na(waiccur) )
{
mlikglob<- mlikcur
mlikcur<- mlikcur
waiccand<- waiccur
waicglob<- waiccur
waiccur<- waiccur
waiccurb<- waiccur
}
}
while(layer>0)
{
withRestarts(tryCatch({
if(printable.opt)print(paste("backward proceed with layer",layer))
if(printable.opt)print(paste("current backward solution is",as.character(varcand)))
vect<-buildmodel(max.cpu = layer,varcur.old = varcurb,statid = model$statid, switch.type = 6,min.N = min.N, max.N = max.N)
if(printable.opt)print(paste("finish backward preparing models at layer",layer))
cluster<-TRUE
flag1<-0
for(mod_id in 1:layer)
{
if(is.null(vect[[mod_id]]$formula))
{
flag1<-flag1+1
}
}
if(flag1==layer)
{
cluster<-FALSE
if(printable.opt)print("!!!!backward Models already estimated!!!!")
}else
{
res.par <- parallelize(X = vect,FUN = .self$fitmodel)
}
if(printable.opt)print(paste("end backward optimizing at layer",layer))
for(mod_id in 1:layer)
{
if(cluster){
fmb<-fm
fm<-res.par[[mod_id]]
waiccand<-Inf
if(is.null(fm)&&(is.na(res.par[[mod_id]]$waic)))
{
varcand<-varcurb
if(printable.opt)print("backward Model Fit Error!?")
next
}
}
varcand<-vect[[mod_id]]$varcur
if(cluster)
{
waiccand<-res.par[[mod_id]]$waic
mlikcand<-res.par[[mod_id]]$mlik
}else if(exists("statistics1"))
{
iidd<-bittodec(varcand)+1
waiccand<-statistics1[iidd,2]
mlikcand<-statistics1[iidd,1]
}else if(exists("hashStat"))
{
iidd<-paste(varcand,collapse = "")
waiccand<-values(hashStat[iidd])[2]
mlikcand<-values(hashStat[iidd])[1]
}
if(objective==0)
{
objcand<-waiccand
objcur<-waiccur
objglob<-waicglob
}else
{
objcand<- -mlikcand
objcur<- -mlikcur
objglob<- -mlikglob
}
if(objcand<=objcur|| mod_id ==1)
{
if(printable.opt)print(paste("backward accept with ", objcand))
objcur<-objcand
varcurb<-varcand
waiccur<-waiccand
varcur<-varcand
mlikcur<-mlikcand
if(objcur<objglob)
{
objglob<-objcur
waicglob<-waiccand
varglob<-varcand
mlikglob<-mlikcand
if(!is.null(fm))
modglob<-fm
if(printable.opt)print(paste("backward global optima with ", objcand))
}
}
}
if(objcur!=objglob)
{
if(model$locstop)
{
break
}else
{
varcurb<-varcur
}
}
}),abort = function(){varcur<-varcurb; fm<-fmb;closeAllConnections();options(error=traceback); onerr<-TRUE})
layer<-layer-1
}
model.prob<-1
model.prob.fix<-1
return(list(varcur = varglob, waiccur = waicglob, mlikcur = mlikglob, log.prob.cur = model.prob,log.prob.fix = model.prob.fix, varglob = varglob, waicglob = waicglob, mlikglob = mlikglob, modglob = modglob))
},
#full selection procedure
full_selection=function(model)
{
if(printable.opt)print(paste("begin full selection procedure!","Careful, ",2^Nvars," models have to be estimated"))
if(Nvars>30)
if(printable.opt)print("Finishing the procedure might well take forever!")
varcand<-array(0,Nvars)
varcurb<-varcand
varglob<-varcand
varglob<-NULL
modglob<-NULL
mlikglob<- model$mlikcur
mlikcur<- model$mlikcur
waiccand<- model$waiccur
waicglob<-model$waiccur
waiccur<-model$waiccur
waiccurb<-model$waiccur
fm<-NULL
fmb<-NULL
ubs<-as.integer(bittodec(array(1,Nvars)) + 1)
ub<-model$ub
totit<-as.integer(ubs/ub) + 1
if(model$totalit<totit)
{
totit<-model$totalit
}
if(is.infinite(mlikcur))
{
vectbg<-buildmodel(max.cpu = 1,varcur.old = varcurb,statid = model$statid,switch.type=8, min.N = min.N,max.N = max.N)
if(!is.null(vectbg[[1]]$formula))
{
bgmod <- lapply(X = vectbg,FUN = .self$fitmodel)
waiccur<-bgmod[[1]]$waic
mlikcur<-bgmod[[1]]$mlik
}
else if(exists("statistics1"))
{
iidd<-bittodec(varcand)+1
waiccur<-statistics1[iidd,2]
mlikcur<-statistics1[iidd,1]
}else if(exists("hashStat"))
{
iidd<-paste(varcand,collapse = "")
waiccur<-values(hashStat[iidd])[2]
mlikcur<-values(hashStat[iidd])[1]
}
if(!is.na(waiccur)&&!is.na(mlikcur))
{
mlikglob<- mlikcur
mlikcand<- mlikcur
waiccand<- waiccur
waicglob<- waiccur
waiccur<- waiccur
waiccurb<- waiccur
}
}
for(i in 1:totit)
{
if(ub*i>ubs)
{
ub<-ubs - ub*(i-1)-1
if(printable.opt)print(paste("last ",ub," iterations to complete"))
varcurb<-varcand
}
withRestarts(tryCatch({
vect<-buildmodel(max.cpu = ub,varcur.old = varcurb,statid = model$statid,switch.type = 7, shift.cpu = model$ub*(i-1),min.N = min.N,max.N = max.N)
if(printable.opt)print(paste("proceed with full ecumeration"))
if(printable.opt)print(paste("current solution is",as.character(varcand)))
cluster<-TRUE
flag1<-0
for(mod_id in 1:ub)
{
if(is.null(vect[[mod_id]]$formula))
{
flag1<-flag1+1
}
}
if(flag1==ub)
{
cluster<-FALSE
if(printable.opt)print("!!!!full models already estimated!!!!")
}else
{
res.par <- parallelize(X = vect,FUN = .self$fitmodel)
}
if(printable.opt)print(paste("end optimizing full ecumeration"))
for(mod_id in 1:ub)
{
if(cluster){
fmb<-fm
fm<-res.par[[mod_id]]
waiccand<-Inf
if(is.null(fm)&&(is.na(res.par[[mod_id]]$waic)))
{
varcand<-varcurb
if(printable.opt)print("full Model Fit Error!?")
next
}
}
varcand<-vect[[mod_id]]$varcur
if(cluster)
{
waiccand<-res.par[[mod_id]]$waic
mlikcand<-res.par[[mod_id]]$mlik
}else if(exists("statistics1"))
{
iidd<-bittodec(varcand)+1
waiccand<-statistics1[iidd,2]
mlikcand<-statistics1[iidd,1]
}else if(exists("hashStat"))
{
iidd<-paste(varcand,collapse = "")
waiccand<-values(hashStat[iidd])[2]
mlikcand<-values(hashStat[iidd])[1]
}
if(objective==0)
{
objcand<-waiccand
objcur<-waiccur
objglob<-waicglob
}else
{
objcand<- -mlikcand
objcur<- -mlikcur
objglob<- -mlikglob
}
if(objcand<=objcur)
{
if(printable.opt)print(paste("full accept with ", objcand))
objcur<-objcand
varcurb<-varcand
waiccur<-waiccand
varcur<-varcand
mlikcur<-mlikcand
if(objcur<objglob)
{
objglob<-objcur
waicglob<-waiccand
varglob<-varcand
mlikglob<-mlikcand
if(!is.null(fm))
modglob<-fm
if(printable.opt)print(paste("full global optima with ", objcand))
}
}
varcurb<-varcand
}
#waiccurb<-waiccur
}),abort = function(){varcur<-varcurb; fm<-fmb;closeAllConnections();options(error=traceback); onerr<-TRUE})
}
# if(length(which(varcur == 1))==(Nvars-1))
# {
# vectbg<-buildmodel(max.cpu = 1,varcur.old = varcurb,statid = model$statid,switch.type=8, min.N = min.N,max.N = max.N)
# if(!is.null(vectbg[[1]]$formula))
# {
# bgmod <- lapply(X = vectbg,FUN = .self$fitmodel)
# waiccand<-bgmod[[1]]$waic
# mlikcand<-bgmod[[1]]$mlik
# }
# else
# {
# iidd<-bittodec(varcur)+1
# waiccand<-statistics1[iidd,2]
# mlikcand<-statistics1[iidd,1]
# }
# if(objective==0)
# {
# objcand<-waiccand
# objcur<-waiccur
# objglob<-waicglob
# }else
# {
# objcand<- -mlikcand
# objcur<- -mlikcur
# objglob<- -mlikglob
# }
#
#
# if(objcand<=objcur)
# {
# if(printable.opt)print(paste("full accept with ", objcand))
# objcur<-objcand
# varcurb<-varcand
# waiccur<-waiccand
# varcur<-varcand
# mlikcur<-mlikcand
#
# if(objcur<objglob)
# {
# objglob<-objcur
# waicglob<-waiccand
# varglob<-varcand
# mlikglob<-mlikcand
# if(!is.null(fm))
# modglob<-fm
#
# if(printable.opt)print(paste("full global optima with ", objcand))
# }
# }
# }
#
model.prob<-1
model.prob.fix<-1
return(list(varcur = varglob, waiccur = waicglob, mlikcur = mlikglob, log.prob.cur = model.prob,log.prob.fix = model.prob.fix, varglob = varglob, waicglob = waicglob, mlikglob = mlikglob, modglob = modglob))
},
#forward backward random dance
forw_backw_walk = function(model)
{
varcur<-rbinom(n = Nvars,size = 1,prob = runif(n = 1,min = 0,max = model$p1))
mlikcur<- -Inf
waiccur<- Inf
for(i in 1:model$steps)
{
fff<-forward_selection(list(varcur=rbinom(n = Nvars,size = 1,prob = runif(n = 1,min = 0,max = model$p1)),mlikcur=-Inf,waiccur =Inf,locstop = FALSE,statid=-1))
if(objective ==1)
{
if(fff$mlikglob>mlikcur)
{
mlikcur<-fff$mlikglob
varcur<-fff$varglob
waiccur<-fff$waicglob
}
}
else
{
if(fff$waicglob<waiccur)
{
mlikcur<-fff$mlikglob
varcur<-fff$varglob
waiccur<-fff$waicglob
}
}
set.seed(i*model$steps)
bbb<-backward_selection(list(varcur=rbinom(n = Nvars,size = 1,prob = runif(n = 1,min = model$p1,max = 1)),mlikcur=-Inf,waiccur =Inf,locstop = FALSE,statid=-1))
if(objective ==1)
{
if(bbb$mlikglob>mlikcur)
{
mlikcur<-bbb$mlikglob
varcur<-bbb$varglob
waiccur<-bbb$waicglob
}
}
else
{
if(bbb$waicglob<waiccur)
{
mlikcur<-bbb$mlikglob
varcur<-bbb$varglob
waiccur<-bbb$waicglob
}
}
}
if(model$reverse == FALSE)
{
vect<-buildmodel(max.cpu = 1,varcur.old = varcur,statid = -1,switch.type = type.randomize,min.N = min.N.randomize,max.N = max.N.randomize)
varcur<-vect[[1]]$varcur
#if(printable.opt)print(varcur)
cluster<-TRUE
if(is.null(vect[[1]]$formula))
{
cluster<-FALSE
if(printable.opt)print("!!!!Back Forw reverse model already estimated!!!!")
}else
{
mod <- lapply(X = vect, FUN = fitmodel)
}
if(cluster)
{
waiccur<-mod[[1]]$waic
mlikcur<-mod[[1]]$mlik
}else if(exists("statistics1"))
{
iidd<-bittodec(varcur)+1
waiccur<-statistics1[iidd,2]
mlikcur<-statistics1[iidd,1]
}else if(exists("hashStat"))
{
iidd<-paste(varcur,collapse = "")
waiccur<-values(hashStat[iidd])[2]
mlikcur<-values(hashStat[iidd])[1]
}
# incorporate what happens for the backward optimization
model.prob<-vect[[1]]$log.mod.switch.prob
model.prob.fix<-vect[[1]]$log.mod.switchback.prob
}else # incorporate what happens for the reverse move
{
model.probs<-calculate.move.logprobabilities(switch.type = type.randomize,varold = varcur, varnew = model$varold,min.N = min.N.randomize,max.N = max.N.randomize)
model.prob<-model.probs$log.switch.forw.prob
model.prob.fix<-model.probs$log.switch.back.prob
}
return(list(varcur = varcur, waiccur = waiccur, mlikcur = mlikcur, log.prob.cur = model.prob,log.prob.fix = model.prob.fix, varglob = varcur, waicglob = waiccur, mlikglob = mlikcur))
},
#local greedy optimization
learnlocalND=function(model)
{
# Step.nd
varcand<-model$varcur
varglob<-model$varcur
varcurb<-model$varcur
mlikcand<- model$mlikcur
waiccand<-model$waiccur
modglob<-NULL
fm<-NULL
fmb<-NULL
opt.achieved<-FALSE
# estimate large jump in a reverse move
if(model$reverse || is.infinite(mlikcand))
{
vectbg<-buildmodel(max.cpu = 1,varcur.old = varcand,statid = model$statid,switch.type=8, min.N = min.N,max.N = max.N)
if(!is.null(vectbg[[1]]$formula))
{
bgmod <- lapply(X = vectbg,FUN = .self$fitmodel)
waiccand<-bgmod[[1]]$waic
mlikcand<-bgmod[[1]]$mlik
}
else if(exists("statistics1"))
{
iidd<-bittodec(varcand)+1
waiccand<-statistics1[iidd,2]
mlikcand<-statistics1[iidd,1]
}else if(exists("hashStat"))
{
iidd<-paste(varcand,collapse = "")
waiccand<-values(hashStat[iidd])[2]
mlikcand<-values(hashStat[iidd])[1]
}
}
if(printable.opt)print(paste("Begin with ",mlikcand))
mlikglob<- mlikcand
mlikcand<- mlikcand
waiccand<- waiccand
waicglob<- waiccand
waiccur<- waiccand
buf.M.nd <- M.nd
if(model$switch.type == 5)
{
buf.M.nd = Nvars - sum(varcurb)
}else if(model$switch.type == 6)
{
buf.M.nd = sum(varcurb)
}
if(buf.M.nd == 0)
buf.M.nd = 1
if(M.nd<buf.M.nd)
buf.M.nd<-M.nd
for(iterat in 1:buf.M.nd)
{
withRestarts(tryCatch({
# statistics <- describe(statistics)
mmax.cpu = max.cpu
if(model$switch.type == 5)
{
mmax.cpu = Nvars - sum(varcurb)
}else if(model$switch.type == 6)
{
mmax.cpu = sum(varcurb)
}
if(mmax.cpu == 0)
mmax.cpu = 1
vect<-buildmodel(max.cpu = mmax.cpu,varcur.old = varcurb,statid = model$statid, switch.type = model$switch.type,min.N = min.N,max.N = max.N)
cluster<-TRUE
flag1<-0
for(mod_id in 1:mmax.cpu)
{
if(is.null(vect[[mod_id]]$formula))
{
flag1<-flag1+1
}
}
if(flag1==mmax.cpu)
{
cluster<-FALSE
if(printable.opt)print("!!!!Greedy Models already estimated!!!!")
}else
{
res.par <- parallelize(X = vect,FUN = .self$fitmodel)
}
for(mod_id in 1:mmax.cpu)
{
varcand1<-vect[[mod_id]]$varcur
if(cluster)
{
waiccand1<-res.par[[mod_id]]$waic
mlikcand1<-res.par[[mod_id]]$mlik
}else if(exists("statistics1"))
{
iidd<-bittodec(varcand1)+1
waiccand1<-statistics1[iidd,2]
mlikcand1<-statistics1[iidd,1]
}else if(exists("hashStat"))
{
iidd<-paste(varcand1,collapse = "")
waiccand1<-values(hashStat[iidd])[2]
mlikcand1<-values(hashStat[iidd])[1]
}
if(objective==0)
{
objcand<-waiccand1
objcur<-waiccand
objglob<-waicglob
}
else
{
objcand<- -mlikcand1
objcur<- -mlikcand
objglob<- -mlikglob
}
if(objcand<objcur || mod_id ==1)
{
varcand<-varcand1
waiccand<-waiccand1
mlikcand<-mlikcand1
if(printable.opt)print(paste("GREEDY update local optima with ", objcand))
if(cluster)
fm<-res.par[[mod_id]]
}
}
#if(printable.opt)print(waiccand)
if(objective==0)
{
objcand<-waiccand1
objcur<-waiccand
objglob<-waicglob
}
else
{
objcand<- -mlikcand1
objcur<- -mlikcand
objglob<- -mlikglob
}
if(objcur<objglob)
{
waicglob<-waiccand
varglob<-varcand
varcurb<-varcand
mlikglob<-mlikcand
if(cluster)
modglob<-fm
if(printable.opt)print(paste("GREEDY update global optima with ", objcur))
}
}),abort = function(){opt.achieved <- TRUE; fm<-fmb;closeAllConnections();options(error=traceback); onerr<-TRUE})
if(objcur!=objglob)
{
if(locstop.nd)
{
break
}else
{
varcurb<-varcand
}
}
}
#!#if(printable.opt)print(points)
if(model$reverse == FALSE)
{
vect<-buildmodel(max.cpu = 1,varcur.old = varcand,statid = model$statid,switch.type = type.randomize,min.N = min.N.randomize,max.N = max.N.randomize)
varcur<-vect[[1]]$varcur
#if(printable.opt)print(varcur)
cluster<-TRUE
if(is.null(vect[[1]]$formula))
{
cluster<-FALSE
if(printable.opt)print("!!!!Greedy reverse model already estimated!!!!")
}else
{
mod <- lapply(X = vect, FUN = fitmodel)
}
if(cluster)
{
waiccur<-mod[[1]]$waic
mlikcur<-mod[[1]]$mlik
}else if(exists("statistics1"))
{
iidd<-bittodec(varcur)+1
waiccur<-statistics1[iidd,2]
mlikcur<-statistics1[iidd,1]
}else if(exists("hashStat"))
{
iidd<-paste(varcur,collapse = "")
waiccur<-values(hashStat[iidd])[2]
mlikcur<-values(hashStat[iidd])[1]
}
# incorporate what happens for the backward optimization
model.prob<-vect[[1]]$log.mod.switch.prob
model.prob.fix<-vect[[1]]$log.mod.switchback.prob
}else # incorporate what happens for the reverse move
{
model.probs<-calculate.move.logprobabilities(switch.type = type.randomize,varold = varcand, varnew = model$varold,min.N = min.N.randomize,max.N = max.N.randomize)
model.prob<-model.probs$log.switch.forw.prob
model.prob.fix<-model.probs$log.switch.back.prob
varcur <- varcand
waiccur <- waiccand
mlikcur <- mlikcand
}
return(list(varcur = varcur, waiccur = waiccur, mlikcur = mlikcur, log.prob.cur = model.prob,log.prob.fix = model.prob.fix, varglob = varglob, waicglob = waicglob, mlikglob = mlikglob, modglob = modglob))
},
#global emjmcmc procedure for model selection (hyper heuristic logic in terms of COP)
modejumping_mcmc=function(glob.model)
{
stm <- proc.time()
if(printable.opt)print("Begin model selection EMJMCMC2016 procedure")
set.seed(runif(n = 1, min = 1, max = seed), kind = NULL, normal.kind = NULL)
acc_moves<-1
accept_old<-1
distrib_of_proposals <- glob.model$distrib_of_proposals
distrib_of_neighbourhoods <- glob.model$distrib_of_neighbourhoods
# do the search and simulations accross the modes
g.results[4,1]<- 0
g.results[4,2]<- 0
forward_selection(list(varcur=rep(0,length(fparam.example)),mlikcur=-Inf,waiccur =Inf,locstop = FALSE,statid=-1))
backward_selection(list(varcur=rep(1,length(fparam.example)),mlikcur=-Inf,waiccur =Inf,locstop = FALSE,statid=-1))
if(exists("statistics1")&&recalc.margin != 2^Nvars)
{
p.add <<- as.array(post_proceed_results(statistics1)$p.post)
}
waiccur<-Inf
waicglob<-Inf
mlikcur<- -Inf
mlikglob<- -Inf
ratcur<- -Inf
#set up initial parameters
if(is.null(glob.model$varcur))
{
if((!is.na(g.results[1,2]))&&g.results[1,2]>0)
{
vec<-dectobit(g.results[1,2]-1)
varcur<-c(array(0,dim = (Nvars -length(vec))),vec)
waiccur<-g.results[2,1]
waicglob<- g.results[2,1]
mlikcur<- g.results[1,1]
mlikglob<- g.results[1,1]
ratcur<- g.results[1,1]
print(paste("initial solution is set with mlik of ",mlikcur))
}else{
vec<-rbinom(n = Nvars,size = 1,prob = 0.5) # generate an initial solution
varcur<-c(array(0,dim = (Nvars -length(vec))),vec)
}
}else if(length(glob.model$varcur[which(glob.model$varcur %in% c(0,1))])==Nvars)
{
varcur<-glob.model$varcur
}
else
{
if(printable.opt)print("Incorrect initial solution set be the user, a random one is generated")
vec<-rbinom(n = Nvars,size = 1,prob = 0.5) # generate an initial solution
varcur<-c(array(0,dim = (Nvars -length(vec))),vec)
}
varcurb<-varcur
varglob<-varcur
modglob<-NULL
p1 = array(data = 0,dim = Nvars)
p2 = array(data = 1,dim = Nvars)
j<-0
p.post<-array(data = 1,dim = Nvars)
waiccur<-Inf
waicglob<-Inf
mlikcur<- -Inf
mlikglob<- -Inf
ratcur<- -Inf
fm<-NULL
eps.emp<-normprob(p1,p2)
max.cpu.buf<-max.cpu
delta.time <- 0
LocImprove<-0
LocNeighbor<-0
max.cpu.buf<-max.cpu.glob
while((eps.emp>=glob.model$eps || j<= glob.model$maxit || j <= glob.model$burnin) && delta.time < glob.model$max.time && g.results[4,1]<= glob.model$trit && g.results[4,2]<= glob.model$trest)
{
p1<-p.post/acc_moves
set.seed(runif(n = 1, min = 1, max = seed*100), kind = NULL, normal.kind = NULL)
LocImprove <- (sample(x = 5,size = 1,prob = distrib_of_proposals) - 1)
LocNeighbor<-(sample(x = 7,size = 1,prob = distrib_of_neighbourhoods[LocImprove+1,]))
switch.type.glob.buf = LocNeighbor
switch.type.buf = LocNeighbor
if(LocNeighbor == 7)
{
switch.type.glob.buf = 9
switch.type.buf = 9
}
#if(printable.opt)print(LocImprove)
j<-j+1
if(j%%glob.model$print.freq == 0)
{
print(paste(j," iterations completed up to now after ",delta.time," cpu minutes"," best MLIK found ",g.results[1,1] ," current mlik found ",mlikcur, "current acceptance ratio ",acc_moves/j))
}
if(j%%100==0)
seed = runif(n = 1,min = 0,max = 100000)
#withRestarts(tryCatch({
varcur<-varcurb
if(LocImprove<=3)
{
vect<-buildmodel(max.cpu = 1,varcur.old = varcurb,statid = 4 + LocImprove,min.N = min.N.glob,max.N = max.N.glob,switch.type = switch.type.glob.buf)
max.cpu.buf = 1
}else
{
vect<-buildmodel(max.cpu = max.cpu.glob,varcur.old = varcurb,statid = 4 + LocImprove,min.N = min.N,max.N = max.N,switch.type = switch.type.glob.buf)
max.cpu.buf = max.cpu.glob
}
cluster<-TRUE
flag1<-0
for(mod_id in 1:max.cpu.buf)
{
if(is.null(vect[[mod_id]]$formula))
{
flag1<-flag1+1
}
}
if(flag1==max.cpu.glob)
{
cluster<-FALSE
if(printable.opt)print("!!!!Models already estimated!!!!")
}else
{
if(max.cpu.glob > 1)
res.par <- parallelize(X = vect,FUN = .self$fitmodel)
else
res.par <- lapply(X = vect,FUN = .self$fitmodel)
}
if(LocImprove>3)
{
if(printable.opt)print("!!!!Proceed with no local improvements!!!!")
p.select.y <- array(data = 0, dim = max.cpu.glob)
for(mod_id in 1:max.cpu.glob)
{
if(cluster)
{
fm<-res.par[[mod_id]]
if(is.null(fm)&&(is.na(res.par[[mod_id]]$waic)))
{
varcand<-varcurb
if(printable.opt)print("GlobMTMCMC Model Fit Error!?")
next
}
}
varcand<-vect[[mod_id]]$varcur
if(cluster)
{
waiccand<-res.par[[mod_id]]$waic
mlikcand<-res.par[[mod_id]]$mlik
}else if(exists("statistics1"))
{
iidd<-bittodec(varcand)+1
waiccand<-statistics1[iidd,2]
mlikcand<-statistics1[iidd,1]
}else if(exists("hashStat"))
{
iidd<-paste(varcand,collapse = "")
waiccand<-values(hashStat[iidd])[2]
mlikcand<-values(hashStat[iidd])[1]
}
if((mlikcand>mlikglob)) #update the parameter of interest
{
if(printable.opt)print(paste("GlobMTMCMC update waic.glob = ", waiccand))
if(printable.opt)print(paste("GlobMTMCMC update waic.glob.mlik = ", mlikglob))
mlikglob<-mlikcand
waicglob<-waiccand
varglob<-varcand
if(cluster)
modglob<-fm
}
g1 <- waiccur
if(waiccur == Inf)
{
g1 = 1
}
p.select.y[mod_id]<-(mlikcand + vect[[mod_id]]$log.mod.switchback.prob+log(lambda(c = cc, alpha = aa, g1 = -g1, g2 = -waiccand,g.domain.pos = FALSE))) # correct for different criteria later
if(is.na(p.select.y[mod_id]))
p.select.y[mod_id] <- 0
if(is.infinite(p.select.y[mod_id]) || p.select.y[mod_id]>100000000)
{
#if(printable.opt)print(paste("very large log.w.y detected ",p.select.y[mod_id]))
p.select.y[mod_id] <- 100000000
}
}
max.p.select.y <- max(p.select.y)
p.select.y<-p.select.y-max.p.select.y
#if(printable.opt)print(paste("max log.w.y is ",max.p.select.y,"normilized log.w.n.y is ", paste(p.select.y,collapse = ", ")))
ID<-sample(x = max.cpu.glob,size = 1,prob = exp(p.select.y))
if(printable.opt)print(paste("cand ",ID," selected"))
varcand<-vect[[ID]]$varcur
if(cluster)
{
waiccand<-res.par[[ID]]$waic
mlikcand<-res.par[[ID]]$mlik
}else if(exists("statistics1"))
{
iidd<-bittodec(varcand)+1
waiccand<-statistics1[iidd,2]
mlikcand<-statistics1[iidd,1]
}else if(exists("hashStat"))
{
iidd<-paste(varcand,collapse = "")
waiccand<-values(hashStat[iidd])[2]
mlikcand<-values(hashStat[iidd])[1]
}
#p.Q.cand<- p.select.y[ID]/sum(p.select.y)
if(printable.opt)print("do reverse step")
p.select.z <- array(data = 0.01, dim = max.cpu.glob)
if(max.cpu.glob!=1)
{
if(switch.type.glob.buf==5)
cstm<-6
else if(switch.type.glob.buf==6)
{
cstm<-5
}
else
{
cstm <- switch.type.glob.buf
}
vect1<-buildmodel(max.cpu = max.cpu.glob -1,varcur.old = varcand,statid = 4 + LocImprove,switch.type = cstm, min.N = min.N, max.N = max.N)
cluster<-TRUE
flag1<-0
for(mod_id in 1:(max.cpu.glob-1))
{
if(is.null(vect1[[mod_id]]$formula))
{
flag1<-flag1+1
}
}
if(flag1==(max.cpu.glob-1))
{
cluster<-FALSE
if(printable.opt)print("!!!!MTMCMC reverse models already estimated!!!!")
}else
{
res.par.back <- parallelize(X = vect1,FUN = .self$fitmodel)
}
for(mod_id in 1:(max.cpu.glob-1))
{
if(cluster)
{
if(is.null(fm)&&(is.na(res.par.back[[mod_id]]$waic)))
{
if(printable.opt)print("locMTMCMC Model Fit Error!?")
next
}
}
varcand.b<-vect1[[mod_id]]$varcur
if(cluster)
{
waiccand.b<-res.par.back[[mod_id]]$waic
mlikcand.b<-res.par.back[[mod_id]]$mlik
}else if(exists("statistics1"))
{
iidd<-bittodec(varcand.b)+1
waiccand.b<-statistics1[iidd,2]
mlikcand.b<-statistics1[iidd,1]
}else if(exists("hashStat"))
{
iidd<-paste(varcand.b,collapse = "")
waiccand.b<-values(hashStat[iidd])[2]
mlikcand.b<-values(hashStat[iidd])[1]
}
if((mlikcand.b>mlikglob))
{
if(printable.opt)print(paste("GlobMTMCMC update waic.glob = ", waiccand.b))
if(printable.opt)print(paste("GlobMTMCMC update waic.glob.mlik = ", mlikcand.b))
mlikglob<-mlikcand.b
waicglob<-waiccand.b
varglob<-varcand.b
if(cluster)
modglob<-fm
}
g1 = waiccand
if(waiccand == Inf)
{
g1 = 1
}
p.select.z[mod_id]<-(mlikcand.b+vect1[[mod_id]]$log.mod.switchback.prob+(lambda(c = cc, alpha = aa, g1 = -g1, g2 = -waiccand.b,g.domain.pos = FALSE))) # correct for different criteria later
if(is.na(p.select.z[mod_id]))
p.select.z[mod_id]=0
if(is.infinite(p.select.z[mod_id]) || p.select.z[mod_id] > 100000000)
{
#if(printable.opt)print(paste("very large log.w.y detected ",p.select.z[mod_id]))
p.select.z[mod_id] <- 100000000
}
}
}
if( waiccur == Inf)
{
g1 = 1
}
p.select.z[max.cpu.glob] <- (mlikcur+vect[[ID]]$log.mod.switch.prob+(lambda(c = cc, alpha = aa, g1 = -g1, g2 = -waiccand,g.domain.pos = FALSE)))
if(is.na(p.select.z[mod_id]))
p.select.z[mod_id]=0
if(is.infinite(p.select.z[mod_id]) || p.select.z[mod_id] > 100000000)
{
#if(printable.opt)print(paste("very large log.w.y detected ",p.select.z[mod_id]))
p.select.z[mod_id] <- 100000000
}
max.p.select.z <- max(p.select.z)
p.select.z<-p.select.z-max.p.select.z
if(printable.opt)print(paste("max log.w.z is ",max.p.select.z,"normilized log.w.n.z is ", paste(p.select.z,collapse = ", ")))
if(log(runif(n = 1,min = 0,max = 1)) < (log(sum(exp(p.select.y)))-log(sum(exp(p.select.z)))) + max.p.select.y - max.p.select.z )
{
mlikcur<-mlikcand
ratcur<-mlikcand
if(printable.opt)print(paste("global MTMCMC update ratcur = ", mlikcand))
if(printable.opt)print(paste("global MTMCMC accept move with ", waiccand))
varcurb<-varcand
waiccur<-waiccand
id<-bittodec(varcurb)+1
acc_moves<-acc_moves+1
if(j<glob.model$burnin)
{
distrib_of_proposals[5]<-distrib_of_proposals[5]+1
}else
{
if(exists("statistics1"))
{
statistics1[id,14]<-statistics1[id,14] + 1
statistics1[id,4]<-statistics1[id,4] + 1
}
p.post<- (p.post + varcurb)
if(allow_offsprings > 0 && acc_moves%%mutation_rate == 0)
{
idmut<-which(varcurb == 1 & p.post > p.allow.tree)
lidmut<-length(idmut)
if(lidmut>0){
imut<-idmut[round(runif(n = 1, min = 1, max = lidmut))]
ltreem<-stri_length(fparam[imut])
if(ltreem<=max.tree.size)
{
crossid<-runif(n = 1,min = 1, max = lidmut)
ltreef<-stri_length(fparam[idmut[crossid]])
if(ltreef+ltreem-4<=max.tree.size)
{
if(allow_offsprings==1)
proposal<-paste(paste(ifelse(runif(n = 1,min = 0,max = 1)<p.nor,"I(","I(1-"),substring(fparam[imut],first=2, last = ltreem),sep = ""),paste(ifelse(runif(n = 1,min = 0,max = 1)<p.nor,"(","(1-"),substring(fparam[idmut[crossid]],first=2, last = ltreef),"))",sep = ""),sep = ifelse(runif(n = 1,min = 0,max = 1)<p.and,"&","|"))
else
proposal<-paste(paste(ifelse(runif(n = 1,min = 0,max = 1)<p.nor,"I(","I(-"),substring(fparam[imut],first=2, last = ltreem),sep = ""),paste("(",substring(fparam[idmut[crossid]],first=2, last = ltreef),"))",sep = ""),sep = ifelse(runif(n = 1,min = 0,max = 1)<p.and,"*","+"))
#maybe check correlations here
if( (!(proposal %in% fparam)) && Nvars<Nvars.max)
{
fparam<<-c(fparam,proposal)
Nvars<<-as.integer(Nvars+1)
print(paste("mutation happended ",proposal," tree added"))
}
else if(!(proposal %in% fparam))
{
to.del<-(which(p.post[(Nvars.init+1):Nvars]<= p.allow.replace) + Nvars.init)
lto.del<-length(x = to.del)
if(lto.del>0)
{
id.replace <- to.del[round(runif(n = 1,min = 1,max = lto.del))]
print(paste("mutation happended ",proposal," tree replaced ", fparam[id.replace]))
fparam[id.replace]<<-proposal
keysarr <- as.array(keys(hashStat))
p.post<-array(data = 0.05,dim = Nvars)
for(jjj in 1:length(keysarr))
{
if(substring(keysarr[jjj],first = id.replace, last = id.replace)=="1")
{
del(x = keysarr[jjj],hash = hashStat)
}
}
}
}
}
}
}
varcurb<-c(varcurb,array(0,dim = (Nvars -length(varcurb))))
varcand<-c(varcand,array(0,dim = (Nvars -length(varcand))))
varglob<-c(varglob,array(0,dim = (Nvars -length(varglob))))
p.post<-c(p.post,array(0,dim = (Nvars -length(p.post))))
p1 = c(p1,array(0,dim = (Nvars -length(p1))))
p2 = c(p1,array(0,dim = (Nvars -length(p1))))
}
}
}else if(j<glob.model$burnin && distrib_of_proposals[5]>0)
{
distrib_of_proposals[5]<-distrib_of_proposals[5] - 1
}
}else
{
if(cluster)
{
fm<-res.par[[mod_id]]
if(is.null(fm)&&(is.na(res.par[[mod_id]]$waic)))
{
varcand<-varcurb
if(printable.opt)print("EMJMCMC Model Fit Error!?")
next
}
}
varcand<-vect[[mod_id]]$varcur
if(cluster)
{
waiccand<-res.par[[mod_id]]$waic
mlikcand<-res.par[[mod_id]]$mlik
}else if(exists("statistics1"))
{
iidd<-bittodec(varcand)+1
waiccand<-statistics1[iidd,2]
mlikcand<-statistics1[iidd,1]
}else if(exists("hashStat"))
{
iidd<-paste(varcand,collapse = "")
waiccand<-values(hashStat[iidd])[2]
mlikcand<-values(hashStat[iidd])[1]
}
varcur<-varcand
}
# try local improvements
if(LocImprove<=3)
{
# sa improvements
if(LocImprove == 0 || LocImprove == 3)
{
if(printable.opt)print("Try SA imptovements")
buf.change <- array(data = 1,dim = Nvars)
buf.change[which(varcur - varcurb!=0)]=0
#buf.opt[floor(runif(n = n.size,min = 1,max = Nvars+0.999999))] = 1
if(objective ==0)
{
objold<-waiccur
objcur<-waiccand
}else
{
objold<- -mlikcur
objcur<- -mlikcand
}
model = list(statid = 4 + LocImprove, switch.type = switch.type.buf, change = buf.change,mlikcur = mlikcur,
varcur = varcur,varold = varcurb, objcur = objcur,objold = objold, sa2 = ifelse(LocImprove == 3,TRUE,FALSE) ,reverse=FALSE)
SA.forw<-learnlocalSA(model)
ratcand<-SA.forw$mlikcur
if(LocImprove == 0)
{
ids = which(buf.change == 1)
model$varcur <- SA.forw$varcur
model$waiccur<-SA.forw$waiccur
model$varcur[ids] = 1 - model$varcur[ids]
model$mlikcur<- -Inf
model$waiccur<- Inf
#estimate the jump!!!
if(objective ==0)
{
model$objold<-waiccur
}else
{
model$objold<- -ratcand
}
model$reverse = TRUE
if(switch.type.buf==5)
{
model$switch.type <- 6
}
else if(switch.type.buf==6)
{
model$switch.type <- 5
}
SA.back<-learnlocalSA(model)
}
if(LocImprove == 0)
{
if(log(runif(n = 1,min = 0,max = 1))<=(ratcand - ratcur - SA.forw$log.prob.cur + SA.forw$log.prob.fix + SA.back$log.prob.cur - SA.back$log.prob.fix))
{
ratcur<-ratcand
mlikcur<-ratcand
if(printable.opt)print(paste("update ratcur through SA = ", ratcur))
varcurb<-SA.forw$varcur
acc_moves<-acc_moves+1
id<-bittodec(varcurb)+1
if(j<glob.model$burnin)
{
distrib_of_proposals[1]<-distrib_of_proposals[1]+1
}else
{
if(exists("statistics1")){
statistics1[id,10]<-statistics1[id,10] + 1
statistics1[id,4]<-statistics1[id,4] + 1
}
p.post<- (p.post + varcurb)
}
}else if(j<glob.model$burnin && distrib_of_proposals[1]>0)
{
distrib_of_proposals[1]<-distrib_of_proposals[1] - 1
}
if((SA.back$mlikglob<mlikglob))
{
if(printable.opt)print(paste("update waic.glob sa.back = ", SA.back$waicglob))
if(printable.opt)print(paste("update waic.glob.mlik sa.back= ", SA.back$mlikglob))
waicglob<-SA.back$waicglob
varglob<-SA.back$varglob
modglob<-SA.back$modglob
}
}else
{
if(log(runif(n = 1,min = 0,max = 1))<=(ratcand - ratcur - SA.forw$log.prob.cur + SA.forw$log.prob.fix + vect[[mod_id]]$log.mod.switchback.prob - vect[[mod_id]]$log.mod.switch.prob))
{
ratcur<-ratcand
mlikcur<-ratcand
if(printable.opt)print(paste("update ratcur through SA = ", ratcur))
varcurb<-SA.forw$varcur
id<-bittodec(varcurb)+1
acc_moves<-acc_moves+1
if(j<glob.model$burnin)
{
distrib_of_proposals[4]<-distrib_of_proposals[4]+1
}else
{
if(exists("statistics1")){
statistics1[id,13]<-statistics1[id,13] + 1
statistics1[id,4]<-statistics1[id,4] + 1
}
p.post<- (p.post + varcurb)
}
}else if(j<glob.model$burnin && distrib_of_proposals[4]>0)
{
distrib_of_proposals[4]<-distrib_of_proposals[4] - 1
}
}
if((SA.forw$mlikglob<mlikglob))
{
if(printable.opt)print(paste("update waic.glob sa.forw = ", SA.forw$waicglob))
if(printable.opt)print(paste("update waic.glob.mlik sa.forw = ", SA.forw$mlikglob))
waicglob<-SA.forw$waicglob
varglob<-SA.forw$varglob
modglob<-SA.forw$modglob
}
}else if(LocImprove == 1)
{
if(printable.opt)print("Try MTMCMC imptovements")
buf.change <- array(data = 1,dim = Nvars)
buf.change[which(varcur - varcurb!=0)]=0
#buf.opt[floor(runif(n = n.size,min = 1,max = Nvars+0.999999))] = 1
model = list(statid = 4 + LocImprove,reverse = FALSE, change = buf.change,mlikcur = mlikcur,waiccur = waiccur,
varcur = varcur,varold = varcurb)
MTMCMC.forw<-learnlocalMCMC(model)
ids = which(buf.change == 1)
model$varcur <- MTMCMC.forw$varcur
model$varcur[ids] = 1 - model$varcur[ids]
model$mlikcur<- -Inf
model$waiccur<- Inf
model$reverse = TRUE
#if(printable.opt)print("learn reverse local MTMCMC")
MTMCMC.back<-learnlocalMCMC(model)
#if(printable.opt)print("finish reverse local MTMCMC")
MTMCMC.p.forw<-MTMCMC.forw$log.prob.cur
MTMCMC.p.back<-MTMCMC.back$log.prob.fix
ratcand<-MTMCMC.forw$mlikcur
if(log(runif(n = 1,min = 0,max = 1))<=(ratcand - ratcur - MTMCMC.forw$log.prob.cur + MTMCMC.forw$log.prob.fix + MTMCMC.back$log.prob.cur - MTMCMC.back$log.prob.fix))
{
ratcur<-ratcand
mlikcur<-ratcand
#if(printable.opt)print(paste("update ratcur through MTMCMC = ", ratcur))
varcurb<-MTMCMC.forw$varcur
acc_moves<-acc_moves+1
id<-bittodec(varcurb)+1
if(j<glob.model$burnin)
{
distrib_of_proposals[2]<-distrib_of_proposals[2]+1
}else
{
if(exists("statistics1")){
statistics1[id,11]<-statistics1[id,11] + 1
statistics1[id,4]<-statistics1[id,4] + 1
}
p.post<- (p.post + varcurb)
}
}else if(j<glob.model$burnin && distrib_of_proposals[2]>1)
{
distrib_of_proposals[2]<-distrib_of_proposals[2] - 1
}
if((MTMCMC.forw$mlikglob<mlikglob))
{
#if(printable.opt)print(paste("update waic.glob MTMCMC.forw = ", MTMCMC.forw$waicglob))
#if(printable.opt)print(paste("update waic.glob.mlik MTMCMC.forw = ", MTMCMC.forw$mlikglob))
waicglob<-MTMCMC.forw$waicglob
varglob<-MTMCMC.forw$varglob
modglob<-MTMCMC.forw$modglob
}
if((MTMCMC.back$mlikglob<mlikglob))
{
#if(printable.opt)print(paste("update waic.glob MTMCMC.back = ", MTMCMC.back$waicglob))
#if(printable.opt)print(paste("update waic.glob.mlik MTMCMC.back= ", MTMCMC.back$mlikglob))
waicglob<-MTMCMC.back$waicglob
varglob<-MTMCMC.back$varglob
modglob<-MTMCMC.back$modglob
}
}else if(LocImprove == 2)
{
if(printable.opt)print("Try greedy heuristic imptovements")
buf.change <- array(data = 1,dim = Nvars)
buf.change[which(varcur - varcurb!=0)]=0
#buf.opt <- buf.change
#buf.opt[floor(runif(n = n.size,min = 1,max = Nvars+0.999999))] = 1
if(objective ==0)
{
objold<-waiccur
objcur<-waiccand
}else
{
objold<- -mlikcur
objcur<- -mlikcand
}
model = list(statid = 4 + LocImprove, change = buf.change,
varcur = varcur,varold = varcurb, switch.type = switch.type.buf,mlikcur = mlikcur,objcur = objcur, objold = objold, reverse=FALSE)
GREEDY.forw<-learnlocalND(model)
ids = which(buf.change == 1)
model$varcur <- GREEDY.forw$varcur
model$varcur[ids] = 1 - model$varcur[ids]
model$mlikcur<- -Inf
model$waiccur<- Inf
model$reverse = TRUE
ratcand<-GREEDY.forw$mlikcur
if(objective ==0)
{
model$objold<-waiccur
}else
{
model$objold<- -ratcand
}
if(switch.type.buf==5)
{
model$switch.type <- 6
}
else if(switch.type.buf==6)
{
model$switch.type <- 5
}
GREEDY.back<-learnlocalND(model)
GREEDY.p.forw<-GREEDY.forw$log.prob.cur
GREEDY.p.back<-GREEDY.back$log.prob.fix
if(log(runif(n = 1,min = 0,max = 1))<=(ratcand - ratcur - GREEDY.forw$log.prob.cur + GREEDY.forw$log.prob.fix + GREEDY.back$log.prob.cur - GREEDY.back$log.prob.fix))
{
ratcur<-ratcand
mlikcur<-ratcand
if(printable.opt)print(paste("update ratcur through ND = ", ratcur))
varcurb<-GREEDY.forw$varcur
acc_moves<-acc_moves+1
id<-bittodec(varcurb)+1
if(j<glob.model$burnin)
{
distrib_of_proposals[3]<-distrib_of_proposals[3]+1
}else
{
if(exists("statistics1")){
statistics1[id,12]<-statistics1[id,12] + 1
statistics1[id,4]<-statistics1[id,4] + 1
}
p.post<- (p.post + varcurb)
}
}else if(j<glob.model$burnin && distrib_of_proposals[3]>1)
{
distrib_of_proposals[3]<-distrib_of_proposals[3] - 1
}
if((GREEDY.forw$mlikglob<mlikglob))
{
if(printable.opt)print(paste("update waic.glob ND.forw = ", GREEDY.forw$waicglob))
if(printable.opt)print(paste("update waic.glob.mlik ND.forw = ", GREEDY.forw$mlikglob))
waicglob<-GREEDY.forw$waicglob
varglob<-GREEDY.forw$varglob
modglob<-GREEDY.forw$modglob
}
if((GREEDY.back$mlikglob<mlikglob))
{
if(printable.opt)print(paste("update waic.glob ND.back = ", GREEDY.back$waicglob))
if(printable.opt)print(paste("update waic.glob.mlik ND.back= ", GREEDY.back$mlikglob))
waicglob<-GREEDY.back$waicglob
varglob<-GREEDY.back$varglob
modglob<-GREEDY.back$modglob
}
}
}
#}),abort = function(){if(printable.opt)print("error");varcur<-varcurb;closeAllConnections();options(error=traceback); onerr<-TRUE})
if(thin_rate!=-1)
{
if(acc_moves == accept_old && j>glob.model$burnin && j%%as.integer(thin_rate)==0) #carry out smart thinning
{
if(!is.null(varcurb))
{
p.post<- (p.post + varcurb)
thin_count <- thin_count +1
id<-bittodec(varcurb)+1
if(exists("statistics1"))
{
statistics1[id,4]<-statistics1[id,4] + 1
}
}
}
}
accept_old <- acc_moves
p2<-p.post/acc_moves
if(j>glob.model$burnin && recalc.margin == 2^Nvars && sum(p2)!=0)
{
p.add <<- p2
}
eps.emp<-normprob(p1,p2)
etm <- proc.time()
delta.time <- (etm[3] - stm[3])/60.0
}
if(is.null(modglob))
{
vect<-buildmodel(max.cpu = 1,varcur.old = varglob,statid = 3, switch.type = 8)
res.par <- lapply(X = vect, FUN = fitmodel)
modglob<-res.par[[1]]$fm
}
#!#print the results
if(printable.opt)print(paste(j," iterations completed in total"))
#if(printable.opt)print(paste("WAIC.glob = ", waicglob))
etm <- proc.time()
tt<-(etm[3]-stm[3])/60.0
acc_ratio <- acc_moves/j
if(printable.opt)print(paste(j," moves proposed in total, ", acc_moves," of them accepted, acceptance ratio is ",acc_ratio))
if(printable.opt)print(paste("posterior distribution ofproposals is", distrib_of_proposals))
id<-bittodec(varglob)+1
if(exists("statistics1"))
{
bayes.res<-post_proceed_results(statistics1)
m.post<-statistics1[,4]/j
}else if(exists("hashStat"))
{
bayes.res<-post_proceed_results_hash(hashStat)
m.post<-NULL
}
else
{
bayes.res<-NULL
m.post<-NULL
}
return(list(model=modglob, vars=varglob, waic = waicglob, mlik = mlikglob, time=(tt), freqs = distrib_of_proposals, acc_ratio = acc_ratio, p.post = p.post/acc_moves, m.post = m.post, p.post.freq = p.post, eps=eps.emp, bayes.results = bayes.res))
},
#save big.data results, if the latter are available
save_results_csv = function(statistics1, filename)
{
write.big.matrix(x = statistics1,filename = filename)
},
#vizualize the results graphically if the latter are available
visualize_results = function(statistics1, template, mds_size, crit,draw_dist = FALSE)
{
ids = NULL
if(crit$mlik){
ids<-c(ids,which(statistics1[,1]==-100000))
ids<-c(ids,which(statistics1[,1]==-Inf))}
#ids<-c(ids,which(statistics1[,1]>0))
if(crit$waic){
ids<-c(ids,which(statistics1[,2]==Inf))
ids<-c(ids,which(statistics1[,3]==Inf))}
if(crit$dic){
ids<-c(ids,which(is.na(statistics1[,3])))}
if(length(ids)!=0)
{
if(crit$mlik)
mlik.lim<-c(min(statistics1[-ids,1],na.rm = TRUE),max(statistics1[-ids,1],na.rm = TRUE))
if(crit$waic)
waic.lim<-c(min(statistics1[-ids,2],na.rm = TRUE),max(statistics1[-ids,2],na.rm = TRUE))
if(crit$dic)
dic.lim<-c(min(statistics1[-ids,3],na.rm = TRUE),max(statistics1[-ids,3],na.rm = TRUE))
}else
{
if(crit$mlik)
mlik.lim<-c(min(statistics1[,1],na.rm = TRUE),max(statistics1[,1],na.rm = TRUE))
if(crit$waic)
waic.lim<-c(min(statistics1[,2],na.rm = TRUE),max(statistics1[,2],na.rm = TRUE))
if(crit$dic)
dic.lim<-c(min(statistics1[,3],na.rm = TRUE),max(statistics1[,3],na.rm = TRUE))
}
norm<-0.5*sqrt(sum(statistics1[,4],na.rm = TRUE))
if(crit$mlik)
{
if(printable.opt)print(paste("drawing ",workdir,template,"_legend.jpg",sep = ""))
capture.output({withRestarts(tryCatch(capture.output({
jpeg(file=paste(workdir,template,"_legend.jpg",sep = ""))
plot(xlab = "posterior probability of a visit found", ylab="visit type",c(0.1,0.2,0.3,0.4,0.5,0.6,0.7),c(7,7,7,-1,-1,-1,-1),ylim = c(0,9), xlim=c(0,1),pch=19, col = 7,cex= c(2,2,2,0,0,0,0))
points(c(0.1,0.2,0.3,0.4,0.5,0.6,0.7),c(6,6,6,-1,-1,-1,-1),pch=8, col = 5,cex= c(1,2,3,0,0,0,0))
points(c(0.1,0.2,0.3,0.4,0.5,0.6,0.7),c(1,1,1,-1,-1,-1,-1),pch=2, col = 2,cex= c(1,2,3,0,0,0,0))
points(c(0.1,0.2,0.3,0.4,0.5,0.6,0.7),c(2,2,2,-1,-1,-1,-1),pch=3, col = 3,cex= c(1,2,3,0,0,0,0))
points(c(0.1,0.2,0.3,0.4,0.5,0.6,0.7),c(3,3,3,-1,-1,-1,-1),pch=4, col = 4,cex= c(1,2,3,0,0,0,0))
points(c(0.1,0.2,0.3,0.4,0.5,0.6,0.7),c(4,4,4,-1,-1,-1,-1),pch=6, col = 6,cex=c(1,2,3,0,0,0,0))
points(c(0.1,0.2,0.3,0.4,0.5,0.6,0.7),c(5,5,5,-1,-1,-1,-1),pch=1, col = 1,cex= c(1,2,3,0,0,0,0))
points(c(0.1,0.2,0.3,0.4,0.5,0.6,0.7),c(8,8,8,-1,-1,-1,-1),pch=19, col = 2,cex= c(2,2,2,0,0,0,0))
legend(x = 0.35,y = 1.5,legend = "- SA 1 local optimization accepted",bty="n" )
legend(x = 0.35,y = 2.5,legend = "- MTMCMC local optimization accepted",bty="n" )
legend(x = 0.35,y = 3.5,legend = "- GREEDY local optimization accepted",bty="n" )
legend(x = 0.35,y = 4.5,legend = "- SA 2 local optimization accepted",bty="n" )
legend(x = 0.35,y = 5.5,legend = "- NO local optimization accepted",bty="n" )
legend(x = 0.35,y = 6.5,legend = "- Totally accepted",bty="n" )
legend(x = 0.35,y = 7.5,legend = "- Totally explored",bty="n" )
legend(x = 0.35,y = 8.5,legend = "- Bayes formula based posterior",bty="n" )
dev.off()
if(printable.opt)print(paste("drawing ",workdir,template,"_mlik.jpg",sep = ""))
jpeg(file=paste(workdir,template,"_mlik.jpg",sep = ""))
plot(ylim = mlik.lim, xlab = "model_id", ylab="MLIK" , statistics1[,1],pch=19, col = 7,cex= 1*((statistics1[,9]+statistics1[,5]+statistics1[,6]+statistics1[,7]+statistics1[,8])>0))
points(statistics1[,1],pch=8, col = ifelse(statistics1[,4]>0,5,0),cex= ifelse(statistics1[,4]>0,statistics1[,4]/norm+1,0))
points(statistics1[,1],pch=2, col = ifelse(statistics1[,10]>0,2,0),cex= ifelse(statistics1[,10]>0,statistics1[,10]/norm+1,0))
points(statistics1[,1],pch=3, col = ifelse(statistics1[,11]>0,3,0),cex= ifelse(statistics1[,11]>0,statistics1[,11]/norm+1,0))
points(statistics1[,1],pch=4, col = ifelse(statistics1[,12]>0,4,0),cex= ifelse(statistics1[,12]>0,statistics1[,12]/norm+1,0))
points(statistics1[,1],pch=6, col = ifelse(statistics1[,13]>0,6,0),cex= ifelse(statistics1[,13]>0,statistics1[,13]/norm+1,0))
points(statistics1[,1],pch=1, col = ifelse(statistics1[,14]>0,1,0),cex= ifelse(statistics1[,14]>0,statistics1[,14]/norm+1,0))
dev.off()
})), abort = function(){onerr<-TRUE})})
}
if(crit$waic)
{
if(printable.opt)print(paste("drawing ",workdir,template,"_waic.jpg",sep = ""))
capture.output({withRestarts(tryCatch(capture.output({
jpeg(file=paste(workdir,template,"_waic.jpg",sep = ""))
plot(ylim = waic.lim, xlab = "model_id", ylab="WAIC" ,statistics1[,2],pch=19, col = 7,cex= 1*((statistics1[,9]+statistics1[,5]+statistics1[,6]+statistics1[,7]+statistics1[,8])>0))
points(statistics1[,2],pch=8, col = ifelse(statistics1[,4]>0,5,0),cex= ifelse(statistics1[,4]>0,statistics1[,4]/norm+1,0))
points(statistics1[,2],pch=2, col = ifelse(statistics1[,10]>0,2,0),cex= ifelse(statistics1[,10]>0,statistics1[,10]/norm+1,0))
points(statistics1[,2],pch=3, col = ifelse(statistics1[,11]>0,3,0),cex= ifelse(statistics1[,11]>0,statistics1[,11]/norm+1,0))
points(statistics1[,2],pch=4, col = ifelse(statistics1[,12]>0,4,0),cex= ifelse(statistics1[,12]>0,statistics1[,12]/norm+1,0))
points(statistics1[,2],pch=6, col = ifelse(statistics1[,13]>0,6,0),cex= ifelse(statistics1[,13]>0,statistics1[,13]/norm+1,0))
points(statistics1[,2],pch=1, col = ifelse(statistics1[,14]>0,1,0),cex= ifelse(statistics1[,14]>0,statistics1[,14]/norm+1,0))
dev.off()
})), abort = function(){onerr<-TRUE})})
}
if(crit$dic)
{
if(printable.opt)print(paste("drawing ",workdir,template,"_dic.jpg",sep = ""))
capture.output({withRestarts(tryCatch(capture.output({
jpeg(file=paste(workdir,template,"_dic.jpg",sep = ""))
plot(ylim = dic.lim, xlab = "model_id", ylab="DIC" ,statistics1[,3],pch=19, col = 7,cex= 1*((statistics1[,9]+statistics1[,5]+statistics1[,6]+statistics1[,7]+statistics1[,8])>0))
points(statistics1[,3],pch=8, col = ifelse(statistics1[,4]>0,5,0),cex= ifelse(statistics1[,4]>0,statistics1[,4]/norm+1,0))
points(statistics1[,3],pch=2, col = ifelse(statistics1[,10]>0,2,0),cex= ifelse(statistics1[,10]>0,statistics1[,10]/norm+1,0))
points(statistics1[,3],pch=3, col = ifelse(statistics1[,11]>0,3,0),cex= ifelse(statistics1[,11]>0,statistics1[,11]/norm+1,0))
points(statistics1[,3],pch=4, col = ifelse(statistics1[,12]>0,4,0),cex= ifelse(statistics1[,12]>0,statistics1[,12]/norm+1,0))
points(statistics1[,3],pch=6, col = ifelse(statistics1[,13]>0,6,0),cex= ifelse(statistics1[,13]>0,statistics1[,13]/norm+1,0))
points(statistics1[,3],pch=1, col = ifelse(statistics1[,14]>0,1,0),cex= ifelse(statistics1[,14]>0,statistics1[,14]/norm+1,0))
dev.off()
})), abort = function(){onerr<-TRUE})})
}
if(crit$waic && crit$mlik)
{
if(printable.opt)print(paste("drawing ",workdir,template,"_mlik-waic.jpg",sep = ""))
capture.output({withRestarts(tryCatch(capture.output({
jpeg(file=paste(workdir,template,"_mlik-waic.jpg",sep = ""))
plot(ylim = mlik.lim, xlim = waic.lim,xlab = "WAIC", ylab="MLIK" ,statistics1[,2],statistics1[,1], pch=19, col = 7,cex= 1*((statistics1[,9]+statistics1[,5]+statistics1[,6]+statistics1[,7]+statistics1[,8])>0))
points(statistics1[,2],statistics1[,1],pch=8, col = ifelse(statistics1[,4]>0,5,0),cex= ifelse(statistics1[,4]>0,statistics1[,4]/norm+1,0))
points(statistics1[,2],statistics1[,1],pch=2, col = ifelse(statistics1[,10]>0,2,0),cex= ifelse(statistics1[,10]>0,statistics1[,10]/norm+1,0))
points(statistics1[,2],statistics1[,1],pch=3, col = ifelse(statistics1[,11]>0,3,0),cex= ifelse(statistics1[,11]>0,statistics1[,11]/norm+1,0))
points(statistics1[,2],statistics1[,1],pch=4, col = ifelse(statistics1[,12]>0,4,0),cex= ifelse(statistics1[,12]>0,statistics1[,12]/norm+1,0))
points(statistics1[,2],statistics1[,1],pch=6, col = ifelse(statistics1[,13]>0,6,0),cex= ifelse(statistics1[,13]>0,statistics1[,13]/norm+1,0))
points(statistics1[,2],statistics1[,1],pch=1, col = ifelse(statistics1[,14]>0,1,0),cex= ifelse(statistics1[,14]>0,statistics1[,14]/norm+1,0))
dev.off()
})), abort = function(){onerr<-TRUE})})
}
if(crit$dic && crit$mlik)
{
if(printable.opt)print(paste("drawing ",workdir,template,"_mlik-dic.jpg",sep = ""))
capture.output({withRestarts(tryCatch(capture.output({
jpeg(file=paste(workdir,template,"_mlik-dic.jpg",sep = ""))
plot(ylim = mlik.lim, xlim = dic.lim,xlab = "DIC", ylab="MLIK" ,statistics1[,3],statistics1[,1], pch=19, col = 7,cex= 1*((statistics1[,9]+statistics1[,5]+statistics1[,6]+statistics1[,7]+statistics1[,8])>0))
points(statistics1[,3],statistics1[,1],pch=8, col = ifelse(statistics1[,4]>0,5,0),cex= ifelse(statistics1[,4]>0,statistics1[,4]/norm+1,0))
points(statistics1[,3],statistics1[,1],pch=2, col = ifelse(statistics1[,10]>0,2,0),cex= ifelse(statistics1[,10]>0,statistics1[,10]/norm+1,0))
points(statistics1[,3],statistics1[,1],pch=3, col = ifelse(statistics1[,11]>0,3,0),cex= ifelse(statistics1[,11]>0,statistics1[,11]/norm+1,0))
points(statistics1[,3],statistics1[,1],pch=4, col = ifelse(statistics1[,12]>0,4,0),cex= ifelse(statistics1[,12]>0,statistics1[,12]/norm+1,0))
points(statistics1[,3],statistics1[,1],pch=6, col = ifelse(statistics1[,13]>0,6,0),cex= ifelse(statistics1[,13]>0,statistics1[,13]/norm+1,0))
points(statistics1[,3],statistics1[,1],pch=1, col = ifelse(statistics1[,14]>0,1,0),cex= ifelse(statistics1[,14]>0,statistics1[,14]/norm+1,0))
dev.off()
})), abort = function(){onerr<-TRUE})})
}
if(crit$dic && crit$waic)
{
if(printable.opt)print(paste("drawing ",workdir,template,"_waic-dic.jpg",sep = ""))
capture.output({withRestarts(tryCatch(capture.output({
jpeg(file=paste(workdir,template,"_waic-dic.jpg",sep = ""))
plot(ylim = waic.lim, xlim = dic.lim,xlab = "WAIC", ylab="DIC" ,statistics1[,3],statistics1[,2], pch=19, col = 7,cex= 1*((statistics1[,9]+statistics1[,5]+statistics1[,6]+statistics1[,7]+statistics1[,8])>0))
points(statistics1[,3],statistics1[,2],pch=8, col = ifelse(statistics1[,4]>0,5,0),cex= ifelse(statistics1[,4]>0,statistics1[,4]/norm+1,0))
points(statistics1[,3],statistics1[,2],pch=2, col = ifelse(statistics1[,10]>0,2,0),cex= ifelse(statistics1[,10]>0,statistics1[,10]/norm+1,0))
points(statistics1[,3],statistics1[,2],pch=3, col = ifelse(statistics1[,11]>0,3,0),cex= ifelse(statistics1[,11]>0,statistics1[,11]/norm+1,0))
points(statistics1[,3],statistics1[,2],pch=4, col = ifelse(statistics1[,12]>0,4,0),cex= ifelse(statistics1[,12]>0,statistics1[,12]/norm+1,0))
points(statistics1[,3],statistics1[,2],pch=6, col = ifelse(statistics1[,13]>0,6,0),cex= ifelse(statistics1[,13]>0,statistics1[,13]/norm+1,0))
points(statistics1[,3],statistics1[,2],pch=1, col = ifelse(statistics1[,14]>0,1,0),cex= ifelse(statistics1[,14]>0,statistics1[,14]/norm+1,0))
dev.off()
})), abort = function(){onerr<-TRUE})})
}
norm1<-(sum(statistics1[,4],na.rm = TRUE))
xyz<-which(!is.na(statistics1[,1]))
xyz<-intersect(xyz,which(statistics1[,1]!=-10000))
moddee<-which(statistics1[,1]==max(statistics1[,1],na.rm = TRUE))[1]
zyx<-array(data = NA,dim = 2 ^(Nvars)+1)
nconsum<-sum(exp(-statistics1[moddee,1]+statistics1[xyz,1]),na.rm = TRUE)
if( nconsum > 0)
{
zyx[xyz]<-exp(statistics1[xyz,1]-statistics1[moddee,1])/nconsum
y.post.lim<-c(0,max(zyx[xyz]))
}else{
zyx[xyz]<-statistics1[xyz,3]/norm1
y.post.lim<-c(0,NaN)
}
if(is.nan(y.post.lim[2]))
y.post.lim[2]<-max(statistics1[,3]/norm1,na.rm = TRUE)
if(is.nan(y.post.lim[2]))
y.post.lim[2]<-1
if(crit$mlik)
{
if(printable.opt)print(paste("drawing ",workdir,template,"_Pr(M|D).jpg",sep = ""))
capture.output({withRestarts(tryCatch(capture.output({
jpeg(file=paste(workdir,template,"_Pr(M|D).jpg",sep = ""))
plot(xlab = "model_id", ylab="Pr(M(model_id)|D)",ylim = y.post.lim, statistics1[,4]/norm1,pch=19, col = 7,cex= 3*((statistics1[,9]+statistics1[,5]+statistics1[,6]+statistics1[,7]+statistics1[,8])>0))
points(statistics1[,4]/norm1,pch=8, col = ifelse(statistics1[,4]>0,5,0),cex= ifelse(statistics1[,4]>0,statistics1[,4]/norm+1,0))
points(statistics1[,4]/norm1,pch=2, col = ifelse(statistics1[,10]>0,2,0),cex= ifelse(statistics1[,10]>0,statistics1[,10]/norm+1,0))
points(statistics1[,4]/norm1,pch=3, col = ifelse(statistics1[,11]>0,3,0),cex= ifelse(statistics1[,11]>0,statistics1[,11]/norm+1,0))
points(statistics1[,4]/norm1,pch=4, col = ifelse(statistics1[,12]>0,4,0),cex= ifelse(statistics1[,12]>0,statistics1[,12]/norm+1,0))
points(statistics1[,4]/norm1,pch=6, col = ifelse(statistics1[,13]>0,6,0),cex= ifelse(statistics1[,13]>0,statistics1[,13]/norm+1,0))
points(statistics1[,4]/norm1,pch=1, col = ifelse(statistics1[,14]>0,1,0),cex= ifelse(statistics1[,14]>0,statistics1[,14]/norm+1,0))
points(zyx,pch=20, col = 2,cex= 2)
dev.off()
})), abort = function(){onerr<-TRUE})})
}
if(crit$waic && crit$mlik)
{
if(printable.opt)print(paste("drawing ",workdir,template,"_waic-Pr(M|D).jpg",sep = ""))
capture.output({withRestarts(tryCatch(capture.output({
jpeg(file=paste(workdir,template,"_waic-Pr(M|D).jpg",sep = ""))
plot(xlim = waic.lim, ylim = y.post.lim,xlab = "WAIC(M)", ylab="Pr(M|D)",statistics1[,2],statistics1[,4]/norm1, pch=19, col = 7,cex= 3*((statistics1[,9]+statistics1[,5]+statistics1[,6]+statistics1[,7]+statistics1[,8])>0))
points(statistics1[,2],statistics1[,4]/norm1,pch=8, col = ifelse(statistics1[,4]>0,5,0),cex= ifelse(statistics1[,4]>0,statistics1[,4]/norm+1,0))
points(statistics1[,2],statistics1[,4]/norm1,pch=2, col = ifelse(statistics1[,10]>0,2,0),cex= ifelse(statistics1[,10]>0,statistics1[,10]/norm+1,0))
points(statistics1[,2],statistics1[,4]/norm1,pch=3, col = ifelse(statistics1[,11]>0,3,0),cex= ifelse(statistics1[,11]>0,statistics1[,11]/norm+1,0))
points(statistics1[,2],statistics1[,4]/norm1,pch=4, col = ifelse(statistics1[,12]>0,4,0),cex= ifelse(statistics1[,12]>0,statistics1[,12]/norm+1,0))
points(statistics1[,2],statistics1[,4]/norm1,pch=6, col = ifelse(statistics1[,13]>0,6,0),cex= ifelse(statistics1[,13]>0,statistics1[,13]/norm+1,0))
points(statistics1[,2],statistics1[,4]/norm1,pch=1, col = ifelse(statistics1[,14]>0,1,0),cex= ifelse(statistics1[,14]>0,statistics1[,14]/norm+1,0))
points(statistics1[,2],zyx,pch=20, col = 2,cex= 2)
dev.off()
})), abort = function(){onerr<-TRUE})})
}
if(crit$dic && crit$mlik)
{
if(printable.opt)print(paste("drawing ",workdir,template,"_dic-Pr(M|D).jpg",sep = ""))
capture.output({withRestarts(tryCatch(capture.output({
jpeg(file=paste(workdir,template,"_dic-Pr(M|D).jpg",sep = ""))
plot(xlim = dic.lim, ylim = y.post.lim,xlab = "dic(M)", ylab="Pr(M|D)",statistics1[,3],statistics1[,4]/norm1, pch=19, col = 7,cex= 3*((statistics1[,9]+statistics1[,5]+statistics1[,6]+statistics1[,7]+statistics1[,8])>0))
points(statistics1[,3],statistics1[,4]/norm1,pch=8, col = ifelse(statistics1[,4]>0,5,0),cex= ifelse(statistics1[,4]>0,statistics1[,4]/norm+1,0))
points(statistics1[,3],statistics1[,4]/norm1,pch=2, col = ifelse(statistics1[,10]>0,2,0),cex= ifelse(statistics1[,10]>0,statistics1[,10]/norm+1,0))
points(statistics1[,3],statistics1[,4]/norm1,pch=3, col = ifelse(statistics1[,11]>0,3,0),cex= ifelse(statistics1[,11]>0,statistics1[,11]/norm+1,0))
points(statistics1[,3],statistics1[,4]/norm1,pch=4, col = ifelse(statistics1[,12]>0,4,0),cex= ifelse(statistics1[,12]>0,statistics1[,12]/norm+1,0))
points(statistics1[,3],statistics1[,4]/norm1,pch=6, col = ifelse(statistics1[,13]>0,6,0),cex= ifelse(statistics1[,13]>0,statistics1[,13]/norm+1,0))
points(statistics1[,3],statistics1[,4]/norm1,pch=1, col = ifelse(statistics1[,14]>0,1,0),cex= ifelse(statistics1[,14]>0,statistics1[,14]/norm+1,0))
points(statistics1[,3],zyx,pch=20, col = 2,cex= 2)
dev.off()
})), abort = function(){onerr<-TRUE})})
}
if(crit$mlik)
{
jpeg(file=paste(workdir,template,"_mlik-Pr(M|D).jpg",sep = ""))
capture.output({withRestarts(tryCatch(capture.output({
plot(xlim = mlik.lim,ylim = y.post.lim, xlab = "MLIK", ylab="Pr(M|D)",statistics1[,1],statistics1[,4]/norm1, pch=19, col = 7,cex= 3*((statistics1[,9]+statistics1[,5]+statistics1[,6]+statistics1[,7]+statistics1[,8])>0))
points(statistics1[,1],statistics1[,4]/norm1,pch=8, col = ifelse(statistics1[,4]>0,5,0),cex= ifelse(statistics1[,4]>0,statistics1[,4]/statistics1[,3]*3 +1,0))
points(statistics1[,1],statistics1[,4]/norm1,pch=2, col = ifelse(statistics1[,10]>0,2,0),cex= ifelse(statistics1[,10]>0,statistics1[,10]/statistics1[,4]*3 +1,0))
points(statistics1[,1],statistics1[,4]/norm1,pch=3, col = ifelse(statistics1[,11]>0,3,0),cex= ifelse(statistics1[,11]>0,statistics1[,11]/statistics1[,4]*3 +1,0))
points(statistics1[,1],statistics1[,4]/norm1,pch=4, col = ifelse(statistics1[,12]>0,4,0),cex= ifelse(statistics1[,12]>0,statistics1[,12]/statistics1[,4]*3 +1,0))
points(statistics1[,1],statistics1[,4]/norm1,pch=6, col = ifelse(statistics1[,13]>0,6,0),cex= ifelse(statistics1[,13]>0,statistics1[,13]/statistics1[,4]*3 +1,0))
points(statistics1[,1],statistics1[,4]/norm1,pch=1, col = ifelse(statistics1[,14]>0,1,0),cex= ifelse(statistics1[,14]>0,statistics1[,14]/statistics1[,4]*3 +1,0))
points(statistics1[,1],zyx,pch=20, col = 2,cex= 2)
dev.off()
})), abort = function(){onerr<-TRUE})})
}
if(draw_dist)
{
if(printable.opt)print("Calculating distance matrix, may take a significant amount of time, may also produce errors if your machine does not have enough memory")
capture.output({withRestarts(tryCatch(capture.output({
lldd<-2^(Nvars)+1
moddee<-which(zyx==max(zyx,na.rm = TRUE))
iidr<-which(statistics1[moddee,1]==max(statistics1[moddee,1],na.rm = TRUE))
iidr<-which(statistics1[moddee[iidr],3]==max(statistics1[moddee[iidr],3],na.rm = TRUE))
moddee<-moddee[iidr]
if(length(moddee)>1)
moddee<-moddee[1]
vec<-dectobit(moddee-1)
varcur<-c(array(0,dim = (Nvars -length(vec))),vec)
df = data.frame(varcur)
for(i in 1:(lldd-1))
{
if(i==moddee)
{
next
}else
{
vec<-dectobit(i-1)
varcur<-c(array(0,dim = (Nvars -length(vec))),vec)
df<-cbind(df,varcur)
#colnames(x = df)[i] <- paste("solution ",i)
}
}
df<-t(df)
x<-dist(x = df,method = "binary")
dists<-c(0,x[1:lldd-1])
#length(dists)
#which(dists==0)
})), abort = function(){onerr<-TRUE})})
if(crit$mlik)
{
if(printable.opt)print(paste("drawing ",workdir,template,"_distance-mlik.jpg",sep = ""))
capture.output({withRestarts(tryCatch(capture.output({
jpeg(file=paste(workdir,template,"_distance-MLIK",sep = ""))
plot(xlab = "x = |M* - M| = distance from the main mode", ylab="MLIK(M)",ylim = mlik.lim,y=c(statistics1[moddee,1],statistics1[-moddee,1]),x=dists,pch=19, col = 7,cex= 1*((statistics1[,9]+statistics1[,5]+statistics1[,6]+statistics1[,7]+statistics1[,8])>0))
points(y=c(statistics1[moddee,1],statistics1[-moddee,1]),x=dists,pch=8, col = ifelse(c(statistics1[moddee,4],statistics1[-moddee,4])>0,5,0),cex= ifelse(c(statistics1[moddee,4],statistics1[-moddee,4])>0,c(statistics1[moddee,4],statistics1[-moddee,4])/norm+1,0))
points(y=c(statistics1[moddee,1],statistics1[-moddee,1]),x=dists,pch=2, col = ifelse(c(statistics1[moddee,10],statistics1[-moddee,10])>0,2,0),cex= ifelse(c(statistics1[moddee,10],statistics1[-moddee,10])>0,c(statistics1[moddee,10],statistics1[-moddee,10])/norm+1,0))
points(y=c(statistics1[moddee,1],statistics1[-moddee,1]),x=dists,pch=3, col = ifelse(c(statistics1[moddee,11],statistics1[-moddee,11])>0,3,0),cex= ifelse(c(statistics1[moddee,11],statistics1[-moddee,11])>0,c(statistics1[moddee,11],statistics1[-moddee,11])/norm+1,0))
points(y=c(statistics1[moddee,1],statistics1[-moddee,1]),x=dists,pch=4, col = ifelse(c(statistics1[moddee,12],statistics1[-moddee,12])>0,4,0),cex= ifelse(c(statistics1[moddee,12],statistics1[-moddee,12])>0,c(statistics1[moddee,12],statistics1[-moddee,12])/norm+1,0))
points(y=c(statistics1[moddee,1],statistics1[-moddee,1]),x=dists,pch=6, col = ifelse(c(statistics1[moddee,13],statistics1[-moddee,13])>0,6,0),cex= ifelse(c(statistics1[moddee,13],statistics1[-moddee,13])>0,c(statistics1[moddee,13],statistics1[-moddee,13])/norm+1,0))
points(y=c(statistics1[moddee,1],statistics1[-moddee,1]),x=dists,pch=1, col = ifelse(c(statistics1[moddee,14],statistics1[-moddee,14])>0,1,0),cex= ifelse(c(statistics1[moddee,14],statistics1[-moddee,14])>0,c(statistics1[moddee,14],statistics1[-moddee,14])/norm+1,0))
dev.off()
})), abort = function(){onerr<-TRUE})})
}
if(crit$waic)
{
if(printable.opt)print(paste("drawing ",workdir,template,"_distance-waic.jpg",sep = ""))
capture.output({withRestarts(tryCatch(capture.output({
jpeg(file=paste(workdir,template,"_distance-waic",sep = ""))
plot(xlab = "x = |M* - M| = distance from the main mode", ylab="WAIC(M)",ylim = waic.lim,y=c(statistics1[moddee,2],statistics1[-moddee,2]),x=dists,pch=19, col = 7,cex= 1*((statistics1[,9]+statistics1[,5]+statistics1[,6]+statistics1[,7]+statistics1[,8])>0))
points(y=c(statistics1[moddee,2],statistics1[-moddee,2]),x=dists,pch=8, col = ifelse(c(statistics1[moddee,4],statistics1[-moddee,4])>0,5,0),cex= ifelse(c(statistics1[moddee,4],statistics1[-moddee,4])>0,c(statistics1[moddee,4],statistics1[-moddee,4])/norm+1,0))
points(y=c(statistics1[moddee,2],statistics1[-moddee,2]),x=dists,pch=2, col = ifelse(c(statistics1[moddee,10],statistics1[-moddee,10])>0,2,0),cex= ifelse(c(statistics1[moddee,10],statistics1[-moddee,10])>0,c(statistics1[moddee,10],statistics1[-moddee,10])/norm+1,0))
points(y=c(statistics1[moddee,2],statistics1[-moddee,2]),x=dists,pch=3, col = ifelse(c(statistics1[moddee,11],statistics1[-moddee,11])>0,3,0),cex= ifelse(c(statistics1[moddee,11],statistics1[-moddee,11])>0,c(statistics1[moddee,11],statistics1[-moddee,11])/norm+1,0))
points(y=c(statistics1[moddee,2],statistics1[-moddee,2]),x=dists,pch=4, col = ifelse(c(statistics1[moddee,12],statistics1[-moddee,12])>0,4,0),cex= ifelse(c(statistics1[moddee,12],statistics1[-moddee,12])>0,c(statistics1[moddee,12],statistics1[-moddee,12])/norm+1,0))
points(y=c(statistics1[moddee,2],statistics1[-moddee,2]),x=dists,pch=6, col = ifelse(c(statistics1[moddee,13],statistics1[-moddee,13])>0,6,0),cex= ifelse(c(statistics1[moddee,13],statistics1[-moddee,13])>0,c(statistics1[moddee,13],statistics1[-moddee,13])/norm+1,0))
points(y=c(statistics1[moddee,2],statistics1[-moddee,2]),x=dists,pch=1, col = ifelse(c(statistics1[moddee,14],statistics1[-moddee,14])>0,1,0),cex= ifelse(c(statistics1[moddee,14],statistics1[-moddee,14])>0,c(statistics1[moddee,14],statistics1[-moddee,14])/norm+1,0))
dev.off()
})), abort = function(){onerr<-TRUE})})
}
if(crit$dic)
{
if(printable.opt)print(paste("drawing ",workdir,template,"_distance-dic.jpg",sep = ""))
capture.output({withRestarts(tryCatch(capture.output({
jpeg(file=paste(workdir,template,"_distance-dic",sep = ""))
plot(xlab = "x = |M* - M| = distance from the main mode", ylab="DIC(M)",ylim = dic.lim,y=c(statistics1[moddee,3],statistics1[-moddee,3]),x=dists,pch=19, col = 7,cex= 1*((statistics1[,9]+statistics1[,5]+statistics1[,6]+statistics1[,7]+statistics1[,8])>0))
points(y=c(statistics1[moddee,3],statistics1[-moddee,3]),x=dists,pch=8, col = ifelse(c(statistics1[moddee,4],statistics1[-moddee,4])>0,5,0),cex= ifelse(c(statistics1[moddee,4],statistics1[-moddee,4])>0,c(statistics1[moddee,4],statistics1[-moddee,4])/norm+1,0))
points(y=c(statistics1[moddee,3],statistics1[-moddee,3]),x=dists,pch=2, col = ifelse(c(statistics1[moddee,10],statistics1[-moddee,10])>0,2,0),cex= ifelse(c(statistics1[moddee,10],statistics1[-moddee,10])>0,c(statistics1[moddee,10],statistics1[-moddee,10])/norm+1,0))
points(y=c(statistics1[moddee,3],statistics1[-moddee,3]),x=dists,pch=3, col = ifelse(c(statistics1[moddee,11],statistics1[-moddee,11])>0,3,0),cex= ifelse(c(statistics1[moddee,11],statistics1[-moddee,11])>0,c(statistics1[moddee,11],statistics1[-moddee,11])/norm+1,0))
points(y=c(statistics1[moddee,3],statistics1[-moddee,3]),x=dists,pch=4, col = ifelse(c(statistics1[moddee,12],statistics1[-moddee,12])>0,4,0),cex= ifelse(c(statistics1[moddee,12],statistics1[-moddee,12])>0,c(statistics1[moddee,12],statistics1[-moddee,12])/norm+1,0))
points(y=c(statistics1[moddee,3],statistics1[-moddee,3]),x=dists,pch=6, col = ifelse(c(statistics1[moddee,13],statistics1[-moddee,13])>0,6,0),cex= ifelse(c(statistics1[moddee,13],statistics1[-moddee,13])>0,c(statistics1[moddee,13],statistics1[-moddee,13])/norm+1,0))
points(y=c(statistics1[moddee,3],statistics1[-moddee,3]),x=dists,pch=1, col = ifelse(c(statistics1[moddee,14],statistics1[-moddee,14])>0,1,0),cex= ifelse(c(statistics1[moddee,14],statistics1[-moddee,14])>0,c(statistics1[moddee,14],statistics1[-moddee,14])/norm+1,0))
dev.off()
})), abort = function(){onerr<-TRUE})})
}
if(crit$mlik)
{
if(printable.opt)print(paste("drawing ",workdir,template,"_distance-Pr(M|D).jpg",sep = ""))
capture.output({withRestarts(tryCatch(capture.output({
jpeg(file=paste(workdir,template,"_distance-Pr(M|D).jpg",sep = ""))
plot(xlab = "x = |M* - M| = distance from the main mode",ylim = y.post.lim, ylab="Pr(M|D)",y=c(statistics1[moddee,4],statistics1[-moddee,4])/norm1,x=dists,pch=19, col = 7,cex= 3*((statistics1[,9]+statistics1[,5]+statistics1[,6]+statistics1[,7]+statistics1[,8])>0))
points(y=c(statistics1[moddee,4],statistics1[-moddee,4])/norm1,x=dists,pch=8, col = ifelse(c(statistics1[moddee,4],statistics1[-moddee,4])>0,5,0),cex= ifelse(c(statistics1[moddee,4],statistics1[-moddee,4])>0,c(statistics1[moddee,4],statistics1[-moddee,4])/norm+1,0))
points(y=c(statistics1[moddee,4],statistics1[-moddee,4])/norm1,x=dists,pch=2, col = ifelse(c(statistics1[moddee,10],statistics1[-moddee,10])>0,2,0),cex= ifelse(c(statistics1[moddee,10],statistics1[-moddee,10])>0,c(statistics1[moddee,10],statistics1[-moddee,10])/norm+1,0))
points(y=c(statistics1[moddee,4],statistics1[-moddee,4])/norm1,x=dists,pch=3, col = ifelse(c(statistics1[moddee,11],statistics1[-moddee,11])>0,3,0),cex= ifelse(c(statistics1[moddee,11],statistics1[-moddee,11])>0,c(statistics1[moddee,11],statistics1[-moddee,11])/norm+1,0))
points(y=c(statistics1[moddee,4],statistics1[-moddee,4])/norm1,x=dists,pch=4, col = ifelse(c(statistics1[moddee,12],statistics1[-moddee,12])>0,4,0),cex= ifelse(c(statistics1[moddee,12],statistics1[-moddee,12])>0,c(statistics1[moddee,12],statistics1[-moddee,12])/norm+1,0))
points(y=c(statistics1[moddee,4],statistics1[-moddee,4])/norm1,x=dists,pch=6, col = ifelse(c(statistics1[moddee,13],statistics1[-moddee,13])>0,6,0),cex= ifelse(c(statistics1[moddee,13],statistics1[-moddee,13])>0,c(statistics1[moddee,13],statistics1[-moddee,13])/norm+1,0))
points(y=c(statistics1[moddee,4],statistics1[-moddee,4])/norm1,x=dists,pch=1, col = ifelse(c(statistics1[moddee,14],statistics1[-moddee,14])>0,1,0),cex= ifelse(c(statistics1[moddee,14],statistics1[-moddee,14])>0,c(statistics1[moddee,14],statistics1[-moddee,14])/norm+1,0))
points(y= c(zyx[moddee],zyx[-moddee]),x=dists,pch=20, col = 2,cex= 2)
dev.off()
})), abort = function(){onerr<-TRUE})})
}
}
if(crit$mlik)
{
if(printable.opt)print(paste("drawing ",workdir,template,"_mds-Pr(M|D).jpg",sep = ""))
if(printable.opt)print("Calculating distance matrix, may take a significant amount of time, may also produce errors if your machine does not have enough memory")
if(printable.opt)print("Calculating mds coordinates, may take a significant amount of time, may also produce errors if your machine does not have enough memory")
capture.output({withRestarts(tryCatch(capture.output({
# further address subset of the set of the best solution of cardinality 1024
if(lldd>mds_size)
{
lldd<-mds_size
quant<-(sort(statistics1[,1],decreasing = TRUE)[lldd+1])
indmds<-which(statistics1[,1]>quant)
length(indmds)
}else{
quant<- -Inf
indmds<-1:(lldd)
}
vec<-dectobit(moddee-1)
varcur<-c(array(0,dim = (Nvars -length(vec))),vec)
df = data.frame(varcur)
for(i in 1:(lldd-1))
{
if(i==moddee)
{
next
}else
{
vec<-dectobit(indmds[i]-1)
varcur<-c(array(0,dim = (Nvars -length(vec))),vec)
df<-cbind(df,varcur)
#colnames(x = df)[i] <- paste("solution ",i)
}
}
df<-t(df)
x<-dist(x = df,method = "binary")
dists<-c(0,x[1:lldd-1])
fit.mds <- cmdscale(d = x,eig=FALSE, k=2) # k is the number of dim
#fit.mds # view results
x.mds <- fit.mds[,1]
y.mds <- fit.mds[,2]
jpeg(file=paste(workdir,template,"_mds_map_posteriors.jpg",sep = ""))
plot(x.mds[], y.mds[], xlab="Coordinate 1", ylab="Coordinate 2", main="Metric MDS", type="p",col=5,pch = 8,cex= ifelse(c(statistics1[moddee,4],statistics1[setdiff(indmds, moddee),4])>0,c(statistics1[moddee,4],statistics1[setdiff(indmds, moddee),4])/norm1*50,0))
points(x.mds[], y.mds[], xlab="Coordinate 1", ylab="Coordinate 2", main="Metric MDS", type="p",col=2,pch = 10,cex= c(zyx[moddee],zyx[setdiff(indmds, moddee)])*50,0)
points(x.mds[], y.mds[], xlab="Coordinate 1", ylab="Coordinate 2", main="Metric MDS", type="p",col=2,pch = 10,cex= c(zyx[moddee],zyx[setdiff(indmds, moddee)])*50,0)
points(x.mds[], y.mds[], xlab="Coordinate 1", ylab="Coordinate 2", main="Metric MDS", type="p",col=5,pch = 8,cex= ifelse(c(statistics1[moddee,4],statistics1[setdiff(indmds, moddee),4])>0,c(statistics1[moddee,4],statistics1[setdiff(indmds, moddee),4])/norm1*50,0))
points(x.mds[], y.mds[], xlab="Coordinate 1", ylab="Coordinate 2", main="Metric MDS", type="p",col=5,pch = 8,cex= ifelse(c(statistics1[moddee,4],statistics1[setdiff(indmds, moddee),4])>0,c(statistics1[moddee,4],statistics1[setdiff(indmds, moddee),4])/norm1*50,0))
points(x.mds[], y.mds[], xlab="Coordinate 1", ylab="Coordinate 2", main="Metric MDS", type="p",col=5,pch = 8,cex= ifelse(c(statistics1[moddee,4],statistics1[setdiff(indmds, moddee),4])>0,c(statistics1[moddee,4],statistics1[setdiff(indmds, moddee),4])/norm1*50,0))
points(x.mds[], y.mds[], xlab="Coordinate 1", ylab="Coordinate 2", main="Metric MDS", type="p",col=5,pch = 8,cex= ifelse(c(statistics1[moddee,4],statistics1[setdiff(indmds, moddee),4])>0,c(statistics1[moddee,4],statistics1[setdiff(indmds, moddee),4])/norm1*50,0))
points(x.mds[], y.mds[], xlab="Coordinate 1", ylab="Coordinate 2", main="Metric MDS", type="p",col=5,pch = 8,cex= ifelse(c(statistics1[moddee,4],statistics1[setdiff(indmds, moddee),4])>0,c(statistics1[moddee,4],statistics1[setdiff(indmds, moddee),4])/norm1*50,0))
points(x.mds[], y.mds[], xlab="Coordinate 1", ylab="Coordinate 2", main="Metric MDS", type="p",col=5,pch = 8,cex= ifelse(c(statistics1[moddee,4],statistics1[setdiff(indmds, moddee),4])>0,c(statistics1[moddee,4],statistics1[setdiff(indmds, moddee),4])/norm1*50,0))
points(x.mds[], y.mds[], xlab="Coordinate 1", ylab="Coordinate 2", main="Metric MDS", type="p",col=5,pch = 8,cex= ifelse(c(statistics1[moddee,4],statistics1[setdiff(indmds, moddee),4])>0,c(statistics1[moddee,4],statistics1[setdiff(indmds, moddee),4])/norm1*50,0))
points(x.mds[], y.mds[], xlab="Coordinate 1", ylab="Coordinate 2", main="Metric MDS", type="p",col=5,pch = 8,cex= ifelse(c(statistics1[moddee,4],statistics1[setdiff(indmds, moddee),4])>0,c(statistics1[moddee,4],statistics1[setdiff(indmds, moddee),4])/norm1*50,0))
points(x.mds[], y.mds[], xlab="Coordinate 1", ylab="Coordinate 2", main="Metric MDS", type="p",col=5,pch = 8,cex= ifelse(c(statistics1[moddee,4],statistics1[setdiff(indmds, moddee),4])>0,c(statistics1[moddee,4],statistics1[setdiff(indmds, moddee),4])/norm1*50,0))
points(x.mds[], y.mds[], xlab="Coordinate 1", ylab="Coordinate 2", main="Metric MDS", type="p",col=5,pch = 8,cex= ifelse(c(statistics1[moddee,4],statistics1[setdiff(indmds, moddee),4])>0,c(statistics1[moddee,4],statistics1[setdiff(indmds, moddee),4])/norm1*50,0))
points(x.mds[], y.mds[], xlab="Coordinate 1", ylab="Coordinate 2", main="Metric MDS", type="p",col=5,pch = 8,cex= ifelse(c(statistics1[moddee,4],statistics1[setdiff(indmds, moddee),4])>0,c(statistics1[moddee,4],statistics1[setdiff(indmds, moddee),4])/norm1*50,0))
points(x.mds[], y.mds[], xlab="Coordinate 1", ylab="Coordinate 2", main="Metric MDS", type="p",col=7,pch = 19,cex= 0.4)
points(x.mds[], y.mds[], xlab="Coordinate 1", ylab="Coordinate 2", main="Metric MDS", type="p",col=1,pch = 19,cex= 0.01)
dev.off()
})), abort = function(){onerr<-TRUE})})
}
},
#calculates posterior probabilities based on a current search
post_proceed_results = function(statistics1)
{
xyz<-which(!is.na(statistics1[,1]))
g.results[4,2] <- length(xyz)
xyz<-intersect(xyz,which(statistics1[,1]!=-10000))
moddee<-which(statistics1[,1]==max(statistics1[,1],na.rm = TRUE))[1]
zyx<-array(data = NA,dim = length(statistics1[,1]))
nconsum<-sum(exp(-statistics1[moddee,1]+statistics1[xyz,1]),na.rm = TRUE)
if( nconsum > 0)
{
zyx[xyz]<-exp(statistics1[xyz,1]-statistics1[moddee,1])/nconsum
}else{
nnnorm<-sum(statistics1[xyz,4],na.rm = T)
if(nnnorm==0)
nnnorm <- 1
zyx[xyz]<-statistics1[xyz,4]/nnnorm
}
statistics1[,15]<-zyx
lldd<-2^(Nvars)+1
p.post<-array(data = 0,dim = Nvars)
for(i in xyz)
{
vec<-dectobit(i-1)
varcur<-c(array(0,dim = (Nvars -length(vec))),vec)
p.post <- (p.post + varcur*statistics1[i,15])
}
if(sum(p.post)==0 || sum(p.post)>Nvars)
{
p.post <- array(data = 0.5,dim = Nvars)
}
return(list(p.post = p.post, m.post = zyx, s.mass = sum(exp(statistics1[xyz,1]),na.rm = TRUE)))
},
post_proceed_results_hash = function(hashStat)
{
lHash<-length(hashStat)
mliks <- values(hashStat)[which((1:(lHash * 3)) %%3 == 1)]
xyz<-which(mliks!=-10000)
g.results[4,2] <- lHash
moddee<-which( mliks ==max( mliks ,na.rm = TRUE))[1]
zyx<-array(data = NA,dim = lHash)
nconsum<-sum(exp(- mliks[moddee]+ mliks[xyz]),na.rm = TRUE)
if( nconsum > 0)
{
zyx[xyz]<-exp(mliks[xyz]- mliks[moddee])/nconsum
}else{
diff<-0-mliks[moddee]
mliks<-mliks+diff
nconsum<-sum(exp(- mliks[moddee]+ mliks[xyz]),na.rm = TRUE)
zyx[xyz]<-exp(mliks[xyz]- mliks[moddee])/nconsum
}
keysarr <- as.array(keys(hashStat))
p.post<-array(data = 0,dim = Nvars)
for(i in 1:lHash)
{
if(is.na(zyx[i]))
next
#vec<-dectobit(strtoi(keysarr[i], base = 0L)-1) # we will have to write some function that would handle laaaaargeee integers!
varcur<- as.integer(strsplit(keysarr[i],split = "")[[1]])
p.post <- (p.post + varcur*zyx[i])
}
if(sum(p.post)==0 || sum(p.post)>Nvars)
{
p.post <- array(data = 0.5,dim = Nvars)
}
return(list(p.post = p.post, m.post = zyx, s.mass = sum(exp(mliks),na.rm = TRUE)))
},
calculate_quality_measures = function(vect,n,truth)
{
rmse.pi <- array(data = 0,Nvars)
bias.pi <- array(data = 0,Nvars)
for(i in 1:n)
{
bias.pi <- (bias.pi + (vect[,i]-truth))
rmse.pi <- (rmse.pi + (vect[,i]^2+truth^2 - 2*vect[,i]*truth))
}
bias.pi <- bias.pi/n
rmse.pi <- rmse.pi/n
return(list(bias.pi = bias.pi, rmse.pi = rmse.pi))
}
)
)
options(bigmemory.typecast.warning=FALSE)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.